
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt             #visualisation
import seaborn as sns   #visualisation

import matplotlib.patches as mpatches
import collections
from sklearn.ensemble import IsolationForest
from scipy.stats import kurtosis, skew, norm, probplot
import pylab
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from pandas.plotting import scatter_matrix
import matplotlib.ticker as ticker
from scipy.optimize import curve_fit

from scipy.stats import mannwhitneyu
from sklearn.cluster import KMeans

from sklearn.decomposition import PCA
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import LabelEncoder, MinMaxScaler, StandardScaler, MaxAbsScaler, RobustScaler
from sklearn.cluster import AgglomerativeClustering

from scipy import stats
from sklearn.metrics import ConfusionMatrixDisplay
from sklearn.metrics import precision_score, recall_score

#Libraries
from sklearn.utils.fixes import loguniform
from sklearn.model_selection import RandomizedSearchCV
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import StratifiedKFold

from sklearn.metrics import classification_report, confusion_matrix  

from sklearn import svm
import sklearn.model_selection as model_selection
from sklearn.metrics import accuracy_score
from sklearn.metrics import f1_score
from sklearn.metrics import mean_squared_error
from sklearn.multiclass import OneVsRestClassifier
from sklearn.multiclass import OneVsOneClassifier

from itertools import cycle
from sklearn.metrics import PrecisionRecallDisplay

from sklearn.metrics import precision_recall_curve
from sklearn.metrics import average_precision_score

from sklearn.base import BaseEstimator, TransformerMixin
# Use label_binarize to be multi-label like settings

from sklearn.preprocessing import LabelBinarizer

from sklearn.model_selection import cross_val_score

from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.feature_selection import SelectFromModel
from sklearn import metrics

from sklearn.metrics import RocCurveDisplay
from sklearn.tree import plot_tree
from sklearn.preprocessing import FunctionTransformer
from sklearn.linear_model import LogisticRegression

import scanpy as sc
import tensorflow as tf
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense
from tensorflow import keras
from scikeras.wrappers import KerasClassifier

import warnings
warnings.filterwarnings("ignore")


H_Smart_unfil = pd.read_csv("HCC1806_SmartS_Unfiltered_Data.txt",delimiter="\ ",engine='python',index_col=0)
H_Smart_metadata= pd.read_csv("HCC1806_SmartS_MetaData.tsv", delimiter = "\t", engine = "python", index_col = 0)


H_Smart_metadata.head()


H_Smart_metadata
print("Dataset dimension:", np.shape(H_Smart_metadata))

# Number of cells with the same plate

plates= collections.Counter(H_Smart_metadata['PCR Plate'])
print('\n')
for x in sorted(plates):
    print(f'nr of cells in plate {x}: ', plates[x])
print('\n')

# Checking cells with different Conditio

cond= collections.Counter(H_Smart_metadata['Condition'])
for x in cond:
    print(f'nr of {x} cells: ', cond[x])

print('\n')
h= collections.Counter(H_Smart_metadata['Hours'])
for x in h:
    print(f'nr of cells that are in the cell line for {x} hours: ', h[x])

# Making the dataset readable
    
print('\nA more readable dataframe:')

# change filename into cell_name
HS_metadata= H_Smart_metadata.copy()
HS_metadata.set_index('Cell name', inplace=True)
HS_metadata.drop(['Pos', 'PreprocessingTag', 'ProcessingComments'], axis=1, inplace=True)
HS_metadata

Dataset dimension: (243, 8)


nr of cells in plate 1:  64
nr of cells in plate 2:  63
nr of cells in plate 3:  68
nr of cells in plate 4:  48


nr of Normo cells:  117
nr of Hypo cells:  126


nr of cells that are in the cell line for 24 hours:  243

A more readable dataframe:


label= HS_metadata['Condition'].copy()
label

Cell name
S123    Normo
S26     Normo
S97      Hypo
S104     Hypo
S4       Hypo
        ...  
S210    Normo
S214    Normo
S199     Hypo
S205    Normo
S236    Normo
Name: Condition, Length: 243, dtype: object


H_Smart_unfil.head()


# Let us make it more readable before going deep into it

def name_clean(data):
    
    # Rename columns and rows in H_Smart_Unfil to remove double quotes
    data = data.rename(columns=lambda x: x.strip('"'), index=lambda x: x.strip('"'))

    column_names = data.columns

    # Define a regular expression pattern to match "S" followed by 1 to 3 digits
    pattern = r'(S\d{1,3})'

    # Extract the "SXXX" part from the column names
    extracted_names = column_names.str.extract(pattern, expand=False)

    # Replace the original column names with the extracted names
    data.columns = extracted_names
    return data
  
H_Smart_unfil = name_clean(H_Smart_unfil)
H_Smart_unfil

H_Smart_metadata= name_clean(H_Smart_metadata.T).T
H_Smart_unfil


print('Dimension of the dataset: ', H_Smart_unfil.shape)
n_rows= H_Smart_unfil.shape[0]
n_col= H_Smart_unfil.shape[1]
print('Number of rows:', n_rows, '\nNumber of columns:', n_col)

Dimension of the dataset:  (23396, 243)
Number of rows: 23396 
Number of columns: 243


if len(H_Smart_unfil.dtypes[H_Smart_unfil.dtypes == 'int64']) == n_col:
    print('All entries are integer numbers')

def empty_cols(df):
    nulldf=[]
    for col in df.columns:
        if df[col].sum()==0:
            nulldf.append(col)
    return nulldf
print('Empty columns:', empty_cols(H_Smart_unfil))
print(f'Empty rows: {empty_cols(H_Smart_unfil.T)}')

All entries are integer numbers
Empty columns: []
Empty rows: []


# Number of zeros for each cells
zeros = pd.Series(np.zeros(len(H_Smart_unfil.columns)), index = H_Smart_unfil.columns)

for col in H_Smart_unfil.columns:
    counter = 0
    for i in H_Smart_unfil[col]:
        if i == 0:
            counter += 1
    zeros[col] = counter
    
plt.hist(zeros, bins = 100, color = 'purple')
plt.title('Count of zero entries')
plt.xlabel('Number of zeros')
plt.ylabel('Number of cells')
plt.show()


H_Smart_unfil.describe()


plt.figure(figsize=(40,10))
sns.violinplot(data= H_Smart_unfil, color = 'purple',cut=0)
plt.setp(plt.gca().get_xticklabels(), rotation= 90)
plt.show()


random_cell = np.random.randint(243)
sns.violinplot(H_Smart_unfil.iloc[:,random_cell], color = 'purple')
plt.title(f'Violin Plot of {H_Smart_unfil.columns[random_cell]}');


# Count depth is the total number of counts per cells
plt.figure(figsize = (40, 10))
H_Smart_unfil.sum().plot(kind = 'bar', color = 'darkviolet')
plt.title('Count Depth')
plt.xlabel('Cells')
plt.ylabel('Number of genes')
plt.show()


# UMI counts is the number of different genes in the dataset
UMI = {}
for x in H_Smart_unfil.columns:
    UMI[f'{x}'] = 0
for x in H_Smart_unfil.columns:
    temp = H_Smart_unfil[x][H_Smart_unfil[x] != 0]
    UMI[f'{x}'] = len(temp)
UMI_df = pd.Series(UMI)
plt.figure(figsize = (40, 10))
UMI_df.plot(kind = 'bar', color='navy')
plt.title('UMI counts')
plt.xlabel('Cells')
plt.ylabel('Count of different cells')
plt.gca().legend([]) 
plt.show()


# Fraction of mitocondrial genes is the ratio of the total number of mitocondrial genes over total genes
mt_genes = []
for i in H_Smart_unfil.index.tolist():
    if i[:2] == 'MT':
        mt_genes.append(i)
d = {}
for i in H_Smart_unfil.columns:
    counter = 0
    for j in H_Smart_unfil[i].index.tolist():
        if j in mt_genes:
            counter += H_Smart_unfil.loc[j, i]
    d[i] = counter / H_Smart_unfil[i].sum()

mt = pd.Series(d)
plt.figure(figsize = (40, 10))
mt.plot(kind = 'bar', color='forestgreen')
plt.title('Mitochondrial genes')
plt.xlabel('Cells')
plt.ylabel('Fraction of mitochondrial genes')
plt.gca().legend([]) 
plt.show()


# Let us plot the count per cells
count= H_Smart_unfil.sum()
plt.hist(count, bins= 100, color = 'purple')
plt.xlabel('Total count')
plt.ylabel('Number of cells')
plt.title('Frequency of gene counts')
plt.show()


colN = H_Smart_unfil.shape[1]
df_skew_cells = []
for i in range(colN) :     
    v_df = H_Smart_unfil[H_Smart_unfil.columns[i]]
    df_skew_cells += [skew(v_df)]   
sns.histplot(df_skew_cells,bins=100, legend = False, color = 'teal')
plt.xlabel('Skewness of single cells expression profile')
plt.show()


df_kurt_cells = []
for i in range(colN) :     
    v_df = H_Smart_unfil[H_Smart_unfil.columns[i]]
    df_kurt_cells += [kurtosis(v_df)]   
sns.histplot(df_kurt_cells,bins=100, color = 'purple')
plt.xlabel('Kurtosis of the data')
plt.show()


plt.figure(figsize=(50, 10))
plt.boxplot(H_Smart_unfil);


plt.figure(figsize=(50, 10))
plt.boxplot(H_Smart_unfil.T);


isolation_forest = IsolationForest(random_state=42)
outlier_pred = isolation_forest.fit_predict(H_Smart_unfil)
outl= H_Smart_unfil.iloc[outlier_pred==-1]
outl


dup_cols = H_Smart_unfil.T[H_Smart_unfil.T.duplicated(keep=False)]
if dup_cols.shape[0] > 0:
    print(f"There are {dup_cols.shape[0]} duplicate columns:\n\n{dup_cols}")
else:
    print("No duplicate columns found")

No duplicate columns found


dup_rows= H_Smart_unfil[H_Smart_unfil.duplicated(keep=False)]
print(f"There are {dup_rows.shape[0]} duplicate genes: \n\n{dup_rows}")

There are 89 duplicate genes: 

           S123  S26  S97  S104  S4  S8  S108  S11  S113  S119  ...  S243  \
MMP23A        0    0    0     0   0   0     0    0     0     0  ...     0   
LINC01647     0    0    0     0   0   0     0    0     0     0  ...     0   
LINC01361     0    0    0     0   0   0     0    0     0     0  ...     0   
ITGA10        0    0    0     0   0   0     0    0     0     0  ...     0   
RORC          0    0    0     0   0   0     0    0     0     0  ...     0   
...         ...  ...  ...   ...  ..  ..   ...  ...   ...   ...  ...   ...   
ENPP7         0    0    0     0   0   0     0    0     0     0  ...     0   
DTNA          0    0    0     0   0   0     0    0     0     0  ...     0   
ALPK2         0    0    0     0   0   0     0    0     0     0  ...     0   
OR7G2         0    0    0     0   0   0     0    0     0     0  ...     0   
PLVAP         0    0    0     0   0   0     0    0     0     0  ...     0   

           S193  S198  S232  S204  S210  S214  S199  S205  S236  
MMP23A        0     0     0     0     0     0     0     0     0  
LINC01647     0     0     0     0     0     0     0     0     0  
LINC01361     0     0     0     0     0     0     0     0     0  
ITGA10        0     0     0     0     0     0     0     0     0  
RORC          0     0     0     0     0     0     0     0     0  
...         ...   ...   ...   ...   ...   ...   ...   ...   ...  
ENPP7         0     0     0     0     0     0     0     0     0  
DTNA          0     0     0     0     0     0     0     0     0  
ALPK2         0     0     0     0     0     0     0     0     0  
OR7G2         0     0     0     0     0     0     0     0     0  
PLVAP         0     0     0     0     0     0     0     0     0  

[89 rows x 243 columns]


plt.figure(figsize = (12,8))
plt.title('Correlation of expression profile between cell lines')
sns.heatmap(H_Smart_unfil.iloc[:, :100].corr(), cmap = 'viridis', fmt ='.01f');
print(f"Number of cell lines plotted: {H_Smart_unfil.iloc[:, :100].corr().shape}")

Number of cell lines plotted: (100, 100)


sm = scatter_matrix(H_Smart_unfil.iloc[:,:20], diagonal = 'kde',figsize = (15,12), color = 'purple')


sns.histplot(H_Smart_unfil.iloc[:, :50].corr(), bins=100, palette = 'viridis')
plt.xlabel("Correlation between samples")
plt.ylabel("Frequency")
plt.legend(loc='upper right')
plt.show()

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.


unfil_c= H_Smart_unfil.corr()
non_corr= unfil_c.columns[(unfil_c.abs().mean() < 0.3)]
print(f"The mean correlation value, on average, is {unfil_c.describe().T['mean'].mean()}")
print("Correlation threshold: 0.3")
print(f"There are {len(non_corr)} cell lines with lower correlation than the threshold:\n\n{non_corr.tolist()} ")

The mean correlation value, on average, is 0.7415974004285173
Correlation threshold: 0.3
There are 9 cell lines with lower correlation than the threshold:

['S98', 'S61', 'S62', 'S157', 'S188', 'S92', 'S168', 'S72', 'S78']


H_Smart_unfil[non_corr.tolist()]


fig, axes= plt.subplots(3,3, figsize=(15,10))
for i, x in enumerate(non_corr.tolist()):
    ax= axes[i//3, i%3]
    ax.hist(H_Smart_unfil[x], bins=30, color = 'purple')
    ax.set_title(f"Distribution of {x}")


plt.figure(figsize = (12,8))
plt.title('Correlation of expression between genes')
sns.heatmap(H_Smart_unfil.T.iloc[:, :500].corr(), cmap = 'viridis', fmt ='.01f');
print(f"Number of Genes plotted: {H_Smart_unfil.T.iloc[:, :500].corr().shape}")

Number of Genes plotted: (500, 500)


plt.figure(figsize=(10,7))
plt.hist(label, bins=5, color='purple', alpha=0.7);
plt.xlabel('Condition', fontsize=12)
plt.ylabel("Frequency", fontsize=12)
plt.title("Distribution of Normoxic and Hypoxic Cell lines")
plt.show()


#add a row to the dataset: 'Condition'
H_Smart_unfil_complete = H_Smart_unfil.T.copy()
H_Smart_unfil_complete['Condition'] = label
H_Smart_unfil_complete = H_Smart_unfil_complete.T
H_Smart_unfil_complete


H_unfil_normo = H_Smart_unfil_complete.T[H_Smart_unfil_complete.T['Condition'] == 'Normo'].T
H_unfil_normo.drop('Condition', axis = 'index', inplace = True)
H_unfil_hypo = H_Smart_unfil_complete.T[H_Smart_unfil_complete.T['Condition'] == 'Hypo'].T
H_unfil_hypo.drop('Condition', axis = 'index', inplace = True)


sns.kdeplot(data=H_unfil_normo.mean(), label = 'Normo', clip = [0,500], color = 'purple')
sns.kdeplot(data=H_unfil_hypo.mean(), label = 'Hypo', clip = [0,500], color = 'forestgreen')
plt.legend()
plt.title('Distribution of Normo e Hypo')
plt.show()


fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(20, 10))

# Plot 1
H_unfil_normo.sum().plot(kind='bar', ax=ax1, color='purple', label = 'Normo')
ax1.set_title('Total counts of genes per cell')
ax1.set_xlabel('Cells')
ax1.set_ylabel('Number of genes')
ax1.legend()

# Plot 2
H_unfil_hypo.sum().plot(kind='bar', ax=ax2, color='forestgreen', label = 'Hypo')
ax2.set_title('Total counts of genes per cell')
ax2.set_xlabel('Cells')
ax2.set_ylabel('Number of genes')
ax2.legend()

plt.tight_layout()
plt.show()


# Perform PCA on the unfiltered dataset
pca = PCA(n_components=2)
pca_result = pca.fit_transform(H_Smart_unfil.T)

# Get condition labels for the samples
conditions = H_Smart_metadata["Condition"]

# Create a PCA plot
plt.figure(figsize=(8, 6))
for condition, color in zip(conditions.unique(), ['purple', 'forestgreen']):
    idx = conditions == condition
    plt.scatter(pca_result[idx, 0], pca_result[idx, 1], c=color, label=condition)
plt.legend()
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('PCA Plot of Gene Expression Data')
plt.show()


# PCA model with three components
pca = PCA(n_components=3)
pca_result = pca.fit_transform(H_Smart_unfil.T)

conditions = H_Smart_metadata["Condition"]


fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111, projection='3d')
for condition, color in zip(conditions.unique(), ['purple', 'forestgreen']):
    idx = conditions == condition
    ax.scatter(pca_result[idx, 0], pca_result[idx, 1], c=color, label=condition)

ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
ax.set_title('PCA 3D Scatter Plot')

plt.show()
# Create a DataFrame with the loading scores for each gene
loading_scores = pd.DataFrame(pca.components_.T, columns=['PC1', 'PC2', 'PC3'], index=H_Smart_unfil.index)

# Find the top 10 contributing genes for each principal component
N = 10 
top_contributing_genes_PC1 = loading_scores.nlargest(N, 'PC1').index
top_contributing_genes_PC2 = loading_scores.nlargest(N, 'PC2').index
top_contributing_genes_PC3 = loading_scores.nlargest(N, 'PC3').index
print(f"Top {N} contributing genes for PC1: {list(top_contributing_genes_PC1)}")
print(f"Top {N} contributing genes for PC2: {list(top_contributing_genes_PC2)}")
print(f"Top {N} contributing genes for PC3: {list(top_contributing_genes_PC3)}")

Top 10 contributing genes for PC1: ['FTL', 'ACTB', 'FTH1', 'LDHA', 'GAPDH', 'BEST1', 'KRT19', 'HSPA5', 'MT-CO1', 'ENO1']
Top 10 contributing genes for PC2: ['FTL', 'FTH1', 'BEST1', 'KRT17', 'AKR1C2', 'KRT13', 'KRT6A', 'KRT16', 'SPP1', 'IGFBP3']
Top 10 contributing genes for PC3: ['KRT19', 'FTH1', 'BEST1', 'LDHA', 'KRT17', 'KRT13', 'GAPDH', 'KRT6A', 'MT-CO1', 'ALDOA']


pca_df = pd.DataFrame(pca_result, columns=['PC1', 'PC2', 'PC3'])
pca_df['Condition'] = H_Smart_unfil.columns.map(H_Smart_metadata['Condition'])
sns.pairplot(pca_df, hue='Condition', vars=['PC1', 'PC2', 'PC3'], palette = 'viridis')
plt.show()


Z1 = linkage(H_Smart_unfil.T, method='ward')
plt.figure(figsize=(40, 20))
plt.title("Hierarchical Clustering Dendrogram of Unfiltered data ")
dendrogram(Z1);


label_series = pd.Series(np.zeros(len(H_Smart_unfil.columns)), index = H_Smart_unfil.columns)
for i in label_series.index:
    if label[i] == 'Normo':
        label_series[i] = int(0)
    if label[i] == 'Hypo':
        label_series[i] = int(1)
label_series = label_series.astype(int)


labels_hier = fcluster(Z1,2, criterion = 'maxclust')

labels_hier = labels_hier -1
counter = 0
First_label = []
Second_label = []
for i in range(len(label_series.to_numpy())):
    if labels_hier[i] == 0:
        First_label.append(counter)
        counter += 1
    elif labels_hier[i] == 1:
        Second_label.append(counter)
        counter += 1


counter1 = 0
counter2 = 0
for i in label_series[label_series == 0].index:
    for j in First_label:
        if i == j:
            counter1 += 1
    for w in Second_label:
        if i == w:
            counter2 += 1
counter1, counter2

(0, 0)


counter1 = 0
counter2 = 0
for i in label_series[label_series == 1].index:
    for j in First_label:
        if i == j:
            counter1 += 1
    for w in Second_label:
        if i == w:
            counter2 += 1
counter1, counter2

(0, 0)


adata = sc.AnnData(H_Smart_unfil.T) # Just transforming into a dataset readable by scanpy


# Calculate the total counts per cell
adata.obs['n_counts']  = adata.X.sum(axis=1)

# Calculate skewness and kurtosis
adata_skewness = skew(adata.obs['n_counts'])
adata_kurtosis = kurtosis(adata.obs['n_counts'])


fig, ax = plt.subplots(figsize=(8, 4))
# Fit a Gaussian distribution to the counts
mu, std = norm.fit(adata.obs['n_counts'])

# Plot the histogram of counts and the fitted Gaussian distribution
sns.distplot(adata.obs['n_counts'], bins=50, kde_kws={'label': 'KDE'},color = 'teal')
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, std)
peak = x[p.argmax()]
sns.lineplot(x=x, y=p, label='Gaussian fit', color = 'purple')

# Add skewness and kurtosis
plt.annotate(f'Skewness: {adata_skewness:.2f}', xy=(0.80, 0.2), xycoords='axes fraction', fontsize=10)
plt.annotate(f'Kurtosis: {adata_kurtosis:.2f}', xy=(0.80, 0.15), xycoords='axes fraction', fontsize=10)

plt.xlabel('Total counts per cell')
plt.ylabel('Frequency')
plt.title('Histogram of total counts per cell')

# Add shaded region from peak to maximum count value
plt.axvspan(xmin=peak, xmax=adata.obs['n_counts'].max(), facecolor='limegreen', alpha=0.2, label = "Parts to analyze")
plt.axvspan(xmin=0, xmax=500000, facecolor='limegreen', alpha=0.2)
plt.legend()
# Format so that scientific notation is not used
ax.xaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=False))
ax.ticklabel_format(style='plain')

plt.show()
print("Estimated peak of counts:", peak)

Estimated peak of counts: 2023939.9672439243


fig, ax = plt.subplots(figsize=(8, 4))
sns.distplot(adata.obs['n_counts'][adata.obs['n_counts']<1000000], kde = True, bins = 30, kde_kws={'label': 'KDE'}, color = 'teal')
plt.axvline(x=475000, color='limegreen', linestyle='--', alpha = 0.2)
plt.axvspan(xmin=0, xmax=475000, facecolor='limegreen', alpha=0.2, label = "Part to delete")

plt.xlabel('Total counts per cell')
plt.ylabel('Frequency')
plt.title('Histogram of total counts per cell')

# Format so that scientific notation is not used
ax.xaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=False))
ax.ticklabel_format(style='plain')
plt.legend()

plt.show()
min_cell = 475000


# Select cells with counts greater than 2023940
high_count_cells = adata[adata.obs['n_counts'] > 2023940]

# Create histogram of counts
counts_hist, counts_bins = np.histogram(high_count_cells.obs['n_counts'], bins=100)

# Define function for Gaussian distribution
def gaussian(x, mu, sigma, A):
    return A * np.exp(-(x - mu)**2 / (2 * sigma**2))

# Estimate initial Gaussian parameters from histogram
A_guess = np.max(counts_hist)

# Fit Gaussian to histogram
popt, pcov = curve_fit(gaussian, counts_bins[:-1], counts_hist, p0=[mu, std, A_guess])
# Estimate upper end of counts based on fitted Gaussian
num_std_dev = 1  # Number of standard deviations from the mean
upper_end = popt[0] + num_std_dev * popt[1]
max_count = np.max(high_count_cells.obs['n_counts'])
upper_end = min(upper_end, max_count)

# Create plot of histogram and fitted Gaussian using seaborn
sns.distplot(high_count_cells.obs['n_counts'], bins=100, kde=False, color = 'teal')
sns.lineplot(x=np.linspace(counts_bins[0], counts_bins[-1], 100), y=gaussian(np.linspace(counts_bins[0], counts_bins[-1], 100), *popt), color='purple', label='Gaussian fit')
plt.axvline(x=upper_end, color='limegreen', linestyle='--', alpha = 0.2, label = "Upper end of counts (1 std. dev.)")
plt.axvspan(xmin=upper_end, xmax=adata.obs["n_counts"].max(), facecolor='limegreen', alpha=0.2, label = "Part to delete")

plt.legend()
plt.show()

print("Estimated upper end of counts (2 standard deviations):", upper_end)
max_cell = upper_end

Estimated upper end of counts (2 standard deviations): 3416472.5503895455


adata_copy = adata.copy()
# Filter cells according to identified QC thresholds:
print('Total number of cells: {:d}'.format(adata_copy.n_obs))

sc.pp.filter_cells(adata_copy, min_counts = min_cell)
print('Number of cells after min count filter: {:d}'.format(adata_copy.n_obs))

sc.pp.filter_cells(adata_copy, max_counts = max_cell)
print('Number of cells after max count filter: {:d}'.format(adata_copy.n_obs))

Total number of cells: 243
Number of cells after min count filter: 229
Number of cells after max count filter: 199


# Calculate new skewness and kurtosis
adata_copy_skewness = skew(adata_copy.obs['n_counts'])
adata_copy_kurtosis = kurtosis(adata_copy.obs['n_counts'])
fig, ax = plt.subplots(figsize=(8, 4))
# Fit a Gaussian distribution to the counts
mu, std = norm.fit(adata_copy.obs['n_counts'])

# Plot the histogram of counts and the fitted Gaussian distribution
sns.distplot(adata_copy.obs['n_counts'], bins=50, kde_kws={'label': 'KDE'}, color = 'teal')
xmin, xmax = plt.xlim()
x = np.linspace(xmin, xmax, 100)
p = norm.pdf(x, mu, std)
peak = x[p.argmax()]
sns.lineplot(x=x, y=p, color='purple', label='Gaussian fit')

# Add skewness and kurtosis
plt.annotate(f'Skewness: {adata_copy_skewness:.2f}', xy=(0.80, 0.2), xycoords='axes fraction', fontsize=10)
plt.annotate(f'Kurtosis: {adata_copy_kurtosis:.2f}', xy=(0.80, 0.15), xycoords='axes fraction', fontsize=10)

plt.xlabel('Total counts per cell')
plt.ylabel('Frequency')
plt.title('Histogram of total counts per cell (Filtered)')

plt.legend()
# Format so that scientific notation is not used
ax.xaxis.set_major_formatter(ticker.ScalarFormatter(useMathText=False))
ax.ticklabel_format(style='plain')
plt.show()


filtered_1_data = adata_copy.to_df()


adata.obs['log_counts'] = np.log(adata.obs['n_counts'])
adata_copy.obs['log_counts'] = np.log(adata_copy.obs['n_counts'])

# Calculate skewness and kurtosis
adata_skewness_log = skew(adata.obs['log_counts'])
adata_kurtosis_log = kurtosis(adata.obs['log_counts'])

# Calculate skewness and kurtosis
adata_copy_skewness_log = skew(adata_copy.obs['log_counts'])
adata_copy_kurtosis_log = kurtosis(adata_copy.obs['log_counts'])


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4))

# Fit a Gaussian distribution to the counts
mu, std = norm.fit(adata.obs['log_counts'])
mu2, std2 = norm.fit(adata_copy.obs['log_counts'])

# Plot the histogram of counts and the fitted Gaussian distribution
sns.distplot(adata.obs['log_counts'], bins=50, kde_kws={'label': 'KDE'}, ax=ax1, color = 'teal')
sns.distplot(adata_copy.obs['log_counts'], bins=50, kde_kws={'label': 'KDE'}, ax=ax2, color = 'teal')

xmin, xmax = ax1.get_xlim()
xmin2, xmax2 = ax2.get_xlim()
x1 = np.linspace(xmin, xmax, 100)
x2 = np.linspace(xmin2, xmax2, 100)
p = norm.pdf(x1, mu, std)
p2 = norm.pdf(x2, mu2, std2)
sns.lineplot(x=x1, y=p, color='purple', label='Gaussian fit', ax=ax1)
sns.lineplot(x=x2, y=p2, color='purple', label='Gaussian fit', ax=ax2)

# Annotate skewness and kurtosis on the plots
ax1.annotate(f'Skewness: {adata_skewness_log:.2f}', xy=(0.0225, 0.7), xycoords='axes fraction', fontsize=10)
ax1.annotate(f'Kurtosis: {adata_kurtosis_log:.2f}', xy=(0.0225, 0.65), xycoords='axes fraction', fontsize=10)

ax2.annotate(f'Skewness: {adata_copy_skewness_log:.2f}', xy=(0.0225, 0.7), xycoords='axes fraction', fontsize=10)
ax2.annotate(f'Kurtosis: {adata_copy_kurtosis_log:.2f}', xy=(0.0225, 0.65), xycoords='axes fraction', fontsize=10)


ax1.set_xlabel('log_counts per cell')
ax1.set_ylabel('Frequency')
ax1.set_title('Histogram of log-total counts per cell (Unfiltered)')
ax1.legend()

ax2.set_xlabel('log_counts per cell')
ax2.set_ylabel('Frequency')
ax2.set_title('Histogram of log-total counts per cell (Filtered)')
ax2.legend()
plt.subplots_adjust(wspace=0.4)

plt.show()


# We define adata_qc1 as the adata without outliers (quality control step 1)
adata_qc1 =sc.AnnData(filtered_1_data)

# Calculate the total counts per cell
adata_qc1.obs['n_counts']  = adata_qc1.X.sum(axis=1)

# Let us isolate the one mitocondiral genes
mt_gene_mask = [gene.startswith(('MT-',"MT")) for gene in adata_qc1.var_names] # Mask for mitochondrial genes
adata_qc1.obs['mt_frac'] = adata_qc1.X[:, mt_gene_mask].sum(1)/adata_qc1.obs['n_counts'] # Mitochondrial frac 
adata_qc1.obs['n_genes'] = (adata_qc1.X > 0).sum(1)


adata_2 = sc.AnnData(H_Smart_unfil.T)
adata_2.obs['n_counts']  = adata_2.X.sum(axis=1)
mt_gene_mask = [gene.startswith(('MT-',"MT")) for gene in adata_2.var_names] # Mask for mitochondrial genes
adata_2.obs['mt_frac'] = adata_2.X[:, mt_gene_mask].sum(1)/adata_2.obs['n_counts'] # Mitochondrial ratio 
adata_2.obs['n_genes'] = (adata_2.X > 0).sum(1)


#Unfiltered
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 5))
ax1.set_title("Unfiltered")
ax2.set_title('Filtered')


plt.subplots_adjust(wspace=0.7)
sns.scatterplot(x = adata_2.obs['n_counts'][adata_2.obs['n_counts']<(3.5*1e6)], y = adata_2.obs['n_genes'][adata_2.obs['n_genes']<13000][adata_2.obs['n_genes']>4000], hue = adata_2.obs['mt_frac'][adata_2.obs['mt_frac']<0.05] ,ax =ax1, palette = 'viridis')
sns.scatterplot(x = adata_qc1.obs['n_counts'], y = adata_qc1.obs['n_genes'][adata_qc1.obs['n_genes']>4000], hue = adata_qc1.obs['mt_frac'], ax = ax2, palette = 'viridis')
# Move the legend of ax1 outside the plot
ax1.legend(loc='center left', bbox_to_anchor=(1, 0.5))

# Move the legend of ax2 outside the plot
ax2.legend(loc='center left', bbox_to_anchor=(1, 0.5))

plt.show()


# Let us see only cells with mitocondrial genes
sns.scatterplot(x = adata_qc1.obs['n_counts'], y = adata_qc1.obs['n_genes'], hue = adata_qc1.obs['mt_frac'], palette = 'viridis')
plt.title("Filtered")
plt.show()


sns.scatterplot(x = adata_qc1.obs['n_genes'], y = adata_qc1.obs['mt_frac'], hue = adata_qc1.obs['n_counts'], palette = 'viridis')

plt.axhline(y=0.05, color='purple', linestyle='--', alpha = 0.2, label = "Threshold")

plt.title("Filtered")
plt.show()


# Filtering of the data
adata_qc2 = adata_qc1[adata_qc1.obs['mt_frac']<0.05]
adata_qc2

View of AnnData object with n_obs × n_vars = 197 × 23396
    obs: 'n_counts', 'mt_frac', 'n_genes'


H_Smart_cf = adata_qc2.to_df()
H_Smart_cf = H_Smart_cf.T


diff =  set(label.index) - set(H_Smart_cf.columns) 
diff = pd.Index(list(diff))


label_cf = label.copy()
label_cf.drop(diff.tolist())
H_Smart_cf_complete = H_Smart_cf.T.copy()
H_Smart_cf_complete['Condition'] = label_cf
H_Smart_cf_complete = H_Smart_cf_complete.T


# As last check we see if there are linearly dependent rows and if so we drop it
corr_matrix_cf = H_Smart_cf.corr()
l = []
for i in H_Smart_cf.columns:
    for j in H_Smart_cf.columns:
        if i != j:
            if abs(corr_matrix_cf.loc[i,j]) == 1 :
                l.append((i,j))
l

[]


#We keep only one of the duplicates by discarding the rest
df_duplicates_number = H_Smart_cf[H_Smart_cf.duplicated(keep = 'first')]
df_duplicates_total = H_Smart_cf[H_Smart_cf.duplicated(keep = False)]


print(f'Total number of duplicate genes: {df_duplicates_number.shape[0]}')
print(f'Number of genes which have a copy: {df_duplicates_total.shape[0] - df_duplicates_number.shape[0]}')

Total number of duplicate genes: 463
Number of genes which have a copy: 161


H_Smart_wd = H_Smart_cf[H_Smart_cf.duplicated(keep = 'first') == False]


zeros = pd.Series(np.zeros(len(df_duplicates_total.index)), index = df_duplicates_total.index)

for ind in df_duplicates_total.index.tolist():
    counter = 0
    for j in df_duplicates_total.T[ind]:
        if j == 0:
            counter += 1
    zeros[ind] = counter
    
plt.hist(zeros, bins = 100, color = 'purple')
plt.title('Zero entries in duplicates')
plt.xlabel('Number of zeros')
plt.ylabel('Number of cells')
plt.show()


H_Stat = H_Smart_wd.T.describe().T[['mean', 'std', 'min','max' ]]
H_Stat


scat = sns.scatterplot(x = 'mean', y = 'std', data= H_Stat, color = 'purple', alpha=0.30)
plt.title('Mean vs Standard ')
plt.show()


H_Smart_wd_log = np.log2(1+ H_Smart_wd)
H_Stat_log = H_Smart_wd_log.T.describe().T[['mean', 'std', 'min','max' ]]
plt.style.use("default")
scat = sns.scatterplot(x = 'mean', y = 'std', data= H_Stat_log, color = 'purple', alpha=0.15)
plt.title('Standard deviation vs Mean')
plt.show()


variance_ratio = H_Stat_log['std'] / H_Stat_log['mean']


sns.histplot(H_Stat_log['std'] / H_Stat_log['mean'], color = 'forestgreen')
plt.title('Variance Ratio histogram')
plt.xlabel('Variance Ratio')
plt.ylabel('Count')
plt.axvline(x=0.7, color='purple', linestyle='dashed', linewidth=1.5)
plt.axvline(x=9, color='purple', linestyle= 'dashed', linewidth=1.5)
plt.show()


def CPM (Data):
    Data2 = Data.copy()
    x = np.zeros(Data.shape[1])
    counter2 = 0
    for i in Data.columns:
        counter = 0
        for j in Data[i]:
            counter += j
        x[counter2] = counter
        counter2 += 1
    counter3 = 0
    for i in Data.columns:
        Data2[i] = np.divide(Data2[i],x[counter3])
        counter3 +=1
    cpm_matrix = Data2 * 1e6
    return cpm_matrix
H_CPM_log = CPM(H_Smart_wd_log)

nonzeros = pd.Series(np.zeros(len(H_CPM_log.index)), index = H_CPM_log.index)
for ind in H_CPM_log.index.tolist():
    counter = 0
    for j in H_CPM_log.T[ind]:
        if j == 0:
            counter += 1
    nonzeros[ind] = 197 - counter
nonzeros = pd.DataFrame(nonzeros, columns = ['nonzeros'])


variance_ratio = pd.DataFrame(variance_ratio, columns = ['variance ratio'])
variance_ratio_counts = variance_ratio.join(nonzeros)


sns.scatterplot(data = variance_ratio_counts, x = 'variance ratio', y = 'nonzeros', color = 'forestgreen', alpha = 0.5 )
plt.axvline(x=0.7, color='purple', linestyle='dashed', linewidth=1.5)
plt.axvline(x=9, color='purple', linestyle= 'dashed', linewidth=1.5)
plt.title('Variance Ratio vs Nonzero ')
plt.show()


variance_ratio_counts_high = variance_ratio_counts[variance_ratio_counts['variance ratio'] < 0.7]
variance_ratio_counts_low = variance_ratio_counts[variance_ratio_counts['variance ratio'] > 9] #Here we are going to zoom on this data to look if they can be removed


sns.scatterplot(data = variance_ratio_counts_high, x = 'variance ratio', y= 'nonzeros', color='forestgreen', alpha=0.15)

<Axes: xlabel='variance ratio', ylabel='nonzeros'>


sns.scatterplot(data = variance_ratio_counts_low, x = 'variance ratio', y = 'nonzeros', color = 'forestgreen' )
plt.title('Data to remove variance Ratio vs NonZeros')
plt.show()


H_var_mean_log = H_Stat_log['std'] /  H_Stat_log['mean']
discard1 = H_var_mean_log[H_var_mean_log < 0.7]
discard2 = H_var_mean_log[H_var_mean_log > 9]


H_Smart_gf1 = H_Smart_wd.copy()
for i in list(discard1.index):
    H_Smart_gf1.drop(f'{i}', inplace = True, axis = 'index')
for i in list(discard2.index):
    H_Smart_gf1.drop(f'{i}', inplace = True, axis = 'index')
H_Smart_gf1


H_Stat_log_gf1 = np.log2(1 + H_Smart_gf1).T.describe().T[['mean', 'std', 'min','max' ]]
plt.style.use("default")
plt.figure(figsize=(6,6))
scat = sns.scatterplot(x = 'mean', y = 'std', data= H_Stat_log_gf1, color = 'purple', alpha = 0.5)
plt.title('Standard deviation vs Mean')
plt.show()


sns.histplot(data = H_Stat_log_gf1['std'] / H_Stat_log_gf1['mean'], color = 'forestgreen')
plt.title('Histogram of Variance Ratio after elimination')
plt.xlabel('Variance Ratio')
plt.show()


# Let us compute and plot the sparsity of the resulting dataset
def count_zero(row): 
    l = len(row)
    c = 0
    for i in range(l): 
        if row[i] == 0: 
            c+=1
    return c
    
sparsity = dict()
sparsity['Sparsity'] = []
sparsity['Genes'] = []

for ind in H_Smart_gf1.index: 
    row = H_Smart_gf1.loc[ind, :]
    s = ( count_zero(row) / float(row.size) )
    sparsity['Genes'].append(ind)
    sparsity['Sparsity'].append(s)


plt.figure(figsize=(11,5))
sns.histplot(data= sparsity['Sparsity'], bins = 100, color = 'forestgreen', alpha = 0.3)
plt.axvline(x = 0.98, color = 'purple', linewidth = 2)
plt.title('Sparsity in the data')
plt.show()


sparsity = pd.Series(sparsity['Sparsity'], index = sparsity['Genes'])
sparsity

WASH7P      0.979695
CICP27      0.954315
DDX11L17    0.822335
WASH9P      0.873096
OR4F29      0.928934
              ...   
MT-TM       0.258883
MT-TS1      0.979695
MT-TD       0.736041
MT-TR       0.482234
MAFIP       0.289340
Length: 12808, dtype: float64


diff = sparsity[sparsity > 0.95].index.tolist()


H_Smart_gf2 = H_Smart_gf1.copy()
for i in diff:
    H_Smart_gf2.drop(i, inplace = True, axis = 'index')
H_Smart_gf2


empty_cols(H_Smart_gf2) #no empty columns

[]


def whitneyu_test(dataset1, dataset2): 
  w_test_list = dict()
  for i in range(dataset1.shape[0]):
    wtest,p = mannwhitneyu(dataset1.iloc[i,:], 
              dataset2.iloc[i, :])
    name = dataset1.index[i]
    w_test_list[name] = p
  sorted_w_test_list = sorted(w_test_list.items(), key=lambda x:x[1])
  return sorted_w_test_list

df_new = H_Smart_gf2.T.join(H_Smart_cf_complete.T['Condition'])

df_normo = df_new[df_new['Condition'] == 'Normo'].iloc[:, :-1]
df_hypo = df_new[df_new['Condition'] == 'Hypo'].iloc[:, :-1]
df_normo = df_normo.T
df_hypo = df_hypo.T

df_normo.head()


wtest = whitneyu_test(df_normo, df_hypo)
w_test_df = pd.DataFrame(wtest, columns = ['Gene', 'P-value']).set_index('Gene')

#plot
sns.histplot(w_test_df, bins = 100, alpha = 0.3, color = 'forestgreen')
plt.axvline(x= 0.05, color = 'purple', linewidth = 2)

<matplotlib.lines.Line2D at 0x18d0895ca00>


significant_genes = w_test_df[w_test_df['P-value'] <= 0.05]
top_15 = w_test_df.sort_values(by = 'P-value').iloc[:15, :]
top_15


def scale(data):
    total_genes_per_cell = data.sum(axis = 0)
    median = np.median(total_genes_per_cell)
    scaling_factor = median / total_genes_per_cell
    data_scaled = data.mul(scaling_factor, axis = 1)
    return data_scaled
H_Filtered_Scaled = scale(H_Smart_gf2)


H_Filtered_Scaled_log = H_Filtered_Scaled.applymap(lambda x: np.log(1+x))


plt.figure(figsize=(40,10))
plot=sns.violinplot(data=H_Filtered_Scaled_log, color="purple",cut=0)
plt.setp(plot.get_xticklabels(), rotation=90)
plt.show()


random_cell = np.random.randint(H_Filtered_Scaled_log.shape[1])
plt.figure(figsize = (20,10))
plot = sns.violinplot(data = H_Filtered_Scaled_log.iloc[:,random_cell], color = 'purple', cut = 0, fit = norm)
plot.set_xlabel(H_Filtered_Scaled_log.columns[random_cell])
plot.set_title('Distribution of a cell over its genes')
plt.show()


# Skewness
colN = np.shape(H_Filtered_Scaled_log)[1]
df_skew_cells = []
for i in range(colN) :     
    v_df = H_Filtered_Scaled_log[H_Filtered_Scaled_log.columns[i]]
    df_skew_cells += [skew(v_df)]   
df_skew_cells
sns.histplot(df_skew_cells,bins=100, legend = False, color = 'purple')
plt.xlabel('Skewness')
plt.ylabel('Count')
plt.title('Skewness of cells in scaled log data')
plt.show()


df_kurt_cells = []
for i in range(colN) :     
    v_df = H_Filtered_Scaled_log[H_Filtered_Scaled_log.columns[i]]
    df_kurt_cells += [kurtosis(v_df)]   
df_kurt_cells
sns.histplot(df_kurt_cells,bins=100, color = 'forestgreen')
plt.xlabel('Kurtosis')
plt.ylabel('Count')
plt.title('Kurtosis of cells in scaled log data')
plt.show()


fig, axs = plt.subplots(3, 1, figsize=(40, 30))

# Plot 1: Count Depth
H_Filtered_Scaled_log.sum().plot(kind = 'bar', ax = axs[0], color = 'darkviolet')
axs[0].set_title('Count Depth')
axs[0].set_xlabel('Cells')
axs[0].set_ylabel('Number of genes')

# Plot 2: UMI Counts

UMI = {}
for x in H_Filtered_Scaled_log.columns:
    UMI[f'{x}'] = 0
for x in H_Filtered_Scaled_log.columns:
    temp = H_Filtered_Scaled_log[x][H_Filtered_Scaled_log[x] != 0]
    UMI[f'{x}'] = len(temp)
UMI_df = pd.Series(UMI)
UMI_df.plot(kind='bar', color='navy', ax=axs[1])
axs[1].set_title('UMI counts')
axs[1].set_xlabel('Cells')
axs[1].set_ylabel('Count of different cells')
axs[1].legend([])

# Plot 3: Mitochondrial genes
mt_genes = []
for i in H_Filtered_Scaled_log.index.tolist():
    if i[:2] == 'MT':
        mt_genes.append(i)
d = {}
for i in H_Filtered_Scaled_log.columns:
    counter = 0
    for j in H_Filtered_Scaled_log[i].index.tolist():
        if j in mt_genes:
            counter += H_Filtered_Scaled_log.loc[j, i]
    d[i] = counter / H_Filtered_Scaled_log[i].sum()
mt = pd.Series(d)
mt.plot(kind='bar', color='forestgreen', ax=axs[2])
axs[2].set_title('Mitochondrial genes')
axs[2].set_xlabel('Cells')
axs[2].set_ylabel('Fraction of mitochondrial genes')
axs[2].legend([])


plt.subplots_adjust(hspace=0.5)

plt.show()


Z2 = linkage(H_Filtered_Scaled_log.T, method='ward')
plt.figure(figsize=(40, 20))
plt.title("Hierarchical Clustering Dendrogram of Filtered Scaled log data")
dendrogram(Z2);


H_Filtered_log = H_Smart_gf2.applymap(lambda x: np.log(1 + x))
mean_values = H_Filtered_log.mean()
H_Filtered_log_mean = H_Filtered_log - mean_values


plt.figure(figsize=(40,10))
plot=sns.violinplot(data= H_Filtered_log_mean, palette="viridis",cut=0)
plt.setp(plot.get_xticklabels(), rotation=90)
plt.show()


random_cell = np.random.randint(H_Filtered_log_mean.shape[1])
plt.figure(figsize = (20,10))
plot = sns.violinplot(data = H_Filtered_log_mean.iloc[:,random_cell], color = 'purple', cut = 0, fit = norm)
plot.set_xlabel(H_Filtered_Scaled_log.columns[random_cell])
plot.set_title('Distribution of a cell over its genes')
plt.show()


colN = np.shape(H_Filtered_log_mean)[1]
colN
df_skew_cells = []
for i in range(colN) :     
    v_df = H_Filtered_log_mean[H_Filtered_log_mean.columns[i]]
    df_skew_cells += [skew(v_df)]   
df_skew_cells
sns.histplot(df_skew_cells,bins=100, legend = False, color = 'purple')
plt.xlabel('Skewness')
plt.ylabel('Count')
plt.title('Skewness of cells in log data with centered mean')
plt.show()


df_kurt_cells = []
for i in range(colN) :     
    v_df = H_Filtered_log_mean[H_Filtered_log_mean.columns[i]]
    df_kurt_cells += [kurtosis(v_df)]   
df_kurt_cells
sns.histplot(df_kurt_cells,bins=100, color = 'forestgreen')
plt.xlabel('Kurtosis')
plt.ylabel('Count')
plt.title('Kurtosis of cells in log data with centered mean')
plt.show()


fig, axs = plt.subplots(3, 1, figsize=(40, 30))

# Plot 1: Count Depth
H_Filtered_log_mean.sum().plot(kind = 'bar', ax = axs[0], color = 'darkviolet')
axs[0].set_title('Count Depth')
axs[0].set_xlabel('Cells')
axs[0].set_ylabel('Number of genes')

# Plot 2: UMI Counts

UMI = {}
for x in H_Filtered_log_mean.columns:
    UMI[f'{x}'] = 0
for x in H_Filtered_log_mean.columns:
    temp = H_Filtered_log_mean[x][H_Filtered_log_mean[x] != 0]
    UMI[f'{x}'] = len(temp)
UMI_df = pd.Series(UMI)
UMI_df.plot(kind='bar', color='navy', ax=axs[1])
axs[1].set_title('UMI counts')
axs[1].set_xlabel('Cells')
axs[1].set_ylabel('Count of different cells')
axs[1].legend([])

# Plot 3: Mitochondrial genes
mt_genes = []
for i in H_Filtered_log_mean.index.tolist():
    if i[:2] == 'MT':
        mt_genes.append(i)
d = {}
for i in H_Filtered_log_mean.columns:
    counter = 0
    for j in H_Filtered_log_mean[i].index.tolist():
        if j in mt_genes:
            counter += H_Filtered_log_mean.loc[j, i]
    d[i] = counter / H_Filtered_log_mean[i].sum()
mt = pd.Series(d)
mt.plot(kind='bar', color='forestgreen', ax=axs[2])
axs[2].set_title('Mitochondrial genes')
axs[2].set_xlabel('Cells')
axs[2].set_ylabel('Fraction of mitochondrial genes')
axs[2].legend([])


plt.subplots_adjust(hspace=0.5)

plt.show()


Z3 = linkage(H_Filtered_log_mean.T, method='ward')
plt.figure(figsize=(40, 20))
plt.title("Hierarchical Clustering Dendrogram of Filtered log data with centered mean")
dendrogram(Z3);


log_mean_df = pd.read_csv(r'C:\Users\marti\OneDrive\AI_Lab\log_and_mean_centering.csv', index_col = False)
scaled_log_df = pd.read_csv(r'C:\Users\marti\OneDrive\AI_Lab\scaled_log.csv')
labels = pd.read_csv(r'C:\Users\marti\OneDrive\AI_Lab\label.csv')

log_mean_df.set_index('Unnamed: 0', inplace = True)
scaled_log_df.set_index('Unnamed: 0', inplace = True)
labels.set_index('Cell name', inplace=True)

log_mean_df = log_mean_df.T.join(labels['Condition'])
scaled_log_df = scaled_log_df.T.join(labels['Condition'])

target = pd.DataFrame({'Normo': [0], 'Hypo': [1]}).T
temp = scaled_log_df.set_index('Condition')
scaled_log_df = temp.join(target)
scaled_log_df.set_index([0], inplace= True)

target = pd.DataFrame({'Normo': [0], 'Hypo': [1]}).T
temp = log_mean_df.set_index('Condition')
log_mean_df = temp.join(target)
log_mean_df.set_index([0], inplace= True)


log_mean_df.head()


scaled_log_df.head()


pipe = PCA(n_components=2, random_state=42)

preprocessed_data = pipe.fit_transform(scaled_log_df)


#Plot
pcadf = pd.DataFrame(preprocessed_data,columns=["component_1", "component_2"])

pcadf["true_label"] = scaled_log_df.index

plt.style.use("default")
plt.figure(figsize=(5, 5))
scat = sns.scatterplot(x = 'component_1', y = 'component_2', s=50, data=pcadf, hue="true_label", palette="viridis")
scat.set_title("PCA in 2 dimensions")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)
plt.show()


pipe = PCA(n_components=2, random_state=42)

preprocessed_data = pipe.fit_transform(log_mean_df)


#Plot
pcadf = pd.DataFrame(preprocessed_data,columns=["component_1", "component_2"])

pcadf["true_label"] = log_mean_df.index

plt.figure(figsize=(5, 5))
scat = sns.scatterplot(x = 'component_1', y = 'component_2', s=50, data=pcadf, hue="true_label", palette="viridis")
scat.set_title("PCA in 2 dimensions")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.0)
plt.show()


pipe = PCA(n_components=3, random_state=42)
preprocessed_data = pipe.fit_transform(log_mean_df)

#Plot

fig = plt.figure(figsize=(5, 5))
ax = fig.add_subplot(111, projection='3d')
for condition, color in zip(log_mean_df.index.unique(), ['limegreen', 'teal']):
    idx = log_mean_df.index == condition
    ax.scatter(preprocessed_data[idx, 0], preprocessed_data[idx, 1], preprocessed_data[idx, 2], c=color, label=condition)

ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
ax.set_zlabel('PC3')
ax.set_title('PCA in 3 dimensions')

plt.show()


k_dataset = log_mean_df

preprocessor = PCA(n_components=2, random_state=42)

clusterer = Pipeline([("kmeans", AgglomerativeClustering(n_clusters=2, linkage = 'complete', metric = 'manhattan'))])

pipe = Pipeline([("preprocessor", preprocessor), ("clusterer", clusterer)])
pipe.fit(k_dataset)
preprocessed_data = pipe["preprocessor"].transform(k_dataset)
predicted_labels = 1 - pipe["clusterer"]["kmeans"].labels_
true_labels = k_dataset.index 

#Plot
pcadf = pd.DataFrame(pipe["preprocessor"].transform(k_dataset),columns=["component_1", "component_2"])

pcadf["predicted_cluster"] = predicted_labels
pcadf["true_label"] = true_labels

plt.figure(figsize=(5, 5))
scat = sns.scatterplot(x = 'component_1', y = 'component_2', s=50, data=pcadf, hue="true_label", style="predicted_cluster",palette="viridis")
scat.set_title("Clustering results for cells")
plt.legend(bbox_to_anchor=(1.05, 1), loc=2,  borderaxespad=0.0)
plt.show()


accuracy_score(pcadf['true_label'], pcadf['predicted_cluster'])

0.9390862944162437


x_pca = []
y_pca = []
a = True
for i in range(1,190):
    pca = PCA(n_components=i, svd_solver = 'full')
    pca_result = pca.fit_transform(log_mean_df)
    x_pca.append(i)
    y_pca.append(pca.explained_variance_ratio_.sum())
fig, ax = plt.subplots(1, figsize=(5,5))

data = pd.DataFrame({
    
    'PC dimensions': x_pca,
    '% explained variance': y_pca
})

sns.lineplot(data=data, x='PC dimensions', y='% explained variance', ax=ax, color = '#04395E')
plt.axvline(x=175, color='limegreen', linestyle='--', alpha = 0.2)
plt.axhline(y=0.9508157355962453, color='limegreen', linestyle='--', alpha = 0.2)

ax.set_title('Number of dimension vs explained variance (first 150 dimensions)')

plt.show()


pipe = PCA(n_components=175, random_state=42)

preprocessed_data = pipe.fit_transform(log_mean_df)

#Plot
pcadf = pd.DataFrame(preprocessed_data)

pcadf["true_label"] = log_mean_df.index

pcadf.set_index("true_label", inplace = True)
pcadf


data = log_mean_df
X_train, X_test, y_train, y_test = model_selection.train_test_split(log_mean_df, log_mean_df.index, train_size=0.80, test_size=0.20, random_state=101)

model = RandomForestClassifier(random_state = 101)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)
accuracy_score(y_test, y_pred)

0.95


rs_space={'max_depth':list(np.arange(10, 100, step=10)) + [None],
              'n_estimators': [1, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500 ],
              'max_features':['sqrt', 'auto'],
              'criterion':['gini','entropy'],
              'min_samples_leaf':np.arange(1,4),
              'min_samples_split':np.arange(2, 10, step=2)
          }
kfold = StratifiedKFold(n_splits=3, shuffle=True, random_state=0)
# Define random search
random_search = RandomizedSearchCV(estimator=RandomForestClassifier(random_state = 101), 
                           param_distributions=rs_space, 
                           n_iter=100,
                           refit='accuracy', 
                           n_jobs=-1, 
                           cv=kfold, 
                           verbose=0)
# Fit grid search
rand = random_search.fit(X_train, y_train)
random_estimator, random_parameters = rand.best_estimator_, rand.best_params_


random_estimator

RandomForestClassifier(max_depth=10, min_samples_leaf=2, n_estimators=150,
                       random_state=101)

RandomForestClassifier(max_depth=10, min_samples_leaf=2, n_estimators=150,
                       random_state=101)


print(random_parameters)

{'n_estimators': 150, 'min_samples_split': 2, 'min_samples_leaf': 2, 'max_features': 'sqrt', 'max_depth': 10, 'criterion': 'gini'}


#Save the classifier in case we run again RandomSearch
random_estimator = RandomForestClassifier(max_depth = 10, min_samples_leaf= 2, n_estimators= 150, random_state= 101).fit(X_train,y_train)


print("Accuracy score:", accuracy_score(y_test, random_estimator.predict(X_test)))
print('Mean squared error:', mean_squared_error(y_test, random_estimator.predict(X_test)))

Accuracy score: 0.975
Mean squared error: 0.025


plt.figure(figsize = (12,5))
plt.plot(random_estimator.feature_importances_, color = 'purple')
# Add labels to the graph
plt.xlabel('Features')
plt.ylabel('Feature Importance Score')
# Add title to the graph
plt.title("Visualizing Important Features")
# Visualize the graph
plt.show()


important_features = pd.DataFrame(random_estimator.feature_importances_, index = log_mean_df.columns).rename(columns = {0: 'Important features'})
important_features = important_features[important_features['Important features'] > 0]
important_features


top_15 = important_features.sort_values(by = 'Important features', ascending = False).iloc[:15,:]
top_15


important_features.shape

(1256, 1)


important_genes = log_mean_df[important_features.index]
important_genes


X_train_temp, X_test_temp, y_train_temp, y_test_temp = model_selection.train_test_split(important_genes, important_genes.index, train_size=0.80, test_size=0.20, random_state=42)


X_train, X_test, y_train, y_test = model_selection.train_test_split(log_mean_df, log_mean_df.index, train_size=0.80, test_size=0.20, random_state=42)


rbf = svm.SVC(kernel='rbf', gamma = 0.01, C = 10 ).fit(X_train, y_train)
poly = svm.SVC(kernel='poly', degree=3, C=1).fit(X_train, y_train)
svm_model = svm.SVC().fit(X_train, y_train)
lin = svm.LinearSVC().fit(X_train, y_train) 

#Prediction
rbf_pred = rbf.predict(X_test)
poly_pred = poly.predict(X_test)
svm_pred = svm_model.predict(X_test)
lin_pred = lin.predict(X_test)


poly_accuracy = accuracy_score(y_test, poly_pred)
poly_f1 = f1_score(y_test, poly_pred,)
poly_error = mean_squared_error(y_test, poly_pred)

print('Accuracy (Polynomial Kernel): ',(poly_accuracy))
print('F1 (Polynomial Kernel): ', poly_f1)
print('Mean squared error: ', poly_error)
print('---------')

rbf_accuracy = accuracy_score(y_test, rbf_pred)
rbf_f1 = f1_score(y_test, rbf_pred)
rbf_error = mean_squared_error(y_test, rbf_pred)

print('Accuracy (RBF Kernel): ',rbf_accuracy)
print('F1 (RBF Kernel): ', rbf_f1)
print('Mean squared error: ', rbf_error)
print('---------')

svm_accuracy = accuracy_score(y_test, svm_pred)
svm_f1 = f1_score(y_test, svm_pred)
svm_error = mean_squared_error(y_test, svm_pred)

print('Accuracy (SVC standard): ',svm_accuracy)
print('F1 (SVC standard): ', svm_f1)
print('Mean squared error standard: ', svm_error)
print('---------')


lin_accuracy = accuracy_score(y_test, lin_pred)
lin_f1 = f1_score(y_test, lin_pred)
lin_error = mean_squared_error(y_test, lin_pred)

print('Accuracy (Linear): ', lin_accuracy)
print('F1 (Linear): ', lin_f1)
print('Mean squared error: ', lin_error)

Accuracy (Polynomial Kernel):  0.775
F1 (Polynomial Kernel):  0.8085106382978724
Mean squared error:  0.225
---------
Accuracy (RBF Kernel):  0.475
F1 (RBF Kernel):  0.6440677966101694
Mean squared error:  0.525
---------
Accuracy (SVC standard):  0.975
F1 (SVC standard):  0.972972972972973
Mean squared error standard:  0.025
---------
Accuracy (Linear):  1.0
F1 (Linear):  1.0
Mean squared error:  0.0


classifier = lin
cross_val = cross_val_score(classifier, X_train, y_train, scoring = 'accuracy', cv = 10)
print("Minimum score:", min(cross_val), "\n", 
      "Maximum score:", max(cross_val), "\n", 
      "Average score:", np.mean(cross_val))

Minimum score: 0.9333333333333333 
 Maximum score: 1.0 
 Average score: 0.9808333333333333


X_train_temp, X_test_temp, y_train_temp, y_test_temp = model_selection.train_test_split(important_genes, important_genes.index, train_size=0.80, test_size=0.20, random_state=42)


poly = svm.SVC(kernel='poly', degree=3, C=1).fit(X_train_temp, y_train_temp)
svm_model = svm.SVC().fit(X_train_temp, y_train_temp)
lin = svm.LinearSVC().fit(X_train_temp, y_train_temp) 

poly_pred = poly.predict(X_test_temp)
svm_pred = svm_model.predict(X_test_temp)
lin_pred = lin.predict(X_test_temp)


poly_accuracy = accuracy_score(y_test_temp, poly_pred)
poly_f1 = f1_score(y_test_temp, poly_pred,)
poly_error = mean_squared_error(y_test_temp, poly_pred)

print('Accuracy (Polynomial Kernel): ',(poly_accuracy))
print('F1 (Polynomial Kernel): ', poly_f1)
print('Mean squared error: ', poly_error)
print('---------')

svm_accuracy = accuracy_score(y_test_temp, svm_pred)
svm_f1 = f1_score(y_test_temp, svm_pred)
svm_error = mean_squared_error(y_test_temp, svm_pred)

print('Accuracy (SVC standard): ',svm_accuracy)
print('F1 (SVC standard): ', svm_f1)
print('Mean squared error standard: ', svm_error)
print('---------')

lin_accuracy = accuracy_score(y_test_temp, lin_pred)
lin_f1 = f1_score(y_test_temp, lin_pred)
lin_error = mean_squared_error(y_test_temp, lin_pred)

print('Accuracy (Linear): ', lin_accuracy)
print('F1 (Linear): ', lin_f1)
print('Mean squared error: ', lin_error)

Accuracy (Polynomial Kernel):  0.85
F1 (Polynomial Kernel):  0.8636363636363636
Mean squared error:  0.15
---------
Accuracy (SVC standard):  1.0
F1 (SVC standard):  1.0
Mean squared error standard:  0.0
---------
Accuracy (Linear):  1.0
F1 (Linear):  1.0
Mean squared error:  0.0


X_train, X_test, y_train, y_test = model_selection.train_test_split(pcadf, pcadf.index, train_size=0.80, test_size=0.20, random_state=42)


poly = svm.SVC(kernel='poly', degree=3, C=1).fit(X_train, y_train)
svm_model = svm.SVC().fit(X_train, y_train)
lin = svm.LinearSVC().fit(X_train, y_train) 

#Prediction
poly_pred = poly.predict(X_test)
svm_pred = svm_model.predict(X_test)
lin_pred = lin.predict(X_test)


poly_accuracy = accuracy_score(y_test, poly_pred)
poly_f1 = f1_score(y_test, poly_pred,)
poly_error = mean_squared_error(y_test, poly_pred)

print('Accuracy (Polynomial Kernel): ',(poly_accuracy))
print('F1 (Polynomial Kernel): ', poly_f1)
print('Mean squared error: ', poly_error)
print('---------')

svm_accuracy = accuracy_score(y_test, svm_pred)
svm_f1 = f1_score(y_test, svm_pred)
svm_error = mean_squared_error(y_test, svm_pred)

print('Accuracy (SVC standard): ',svm_accuracy)
print('F1 (SVC standard): ', svm_f1)
print('Mean squared error standard: ', svm_error)
print('---------')


lin_accuracy = accuracy_score(y_test, lin_pred)
lin_f1 = f1_score(y_test, lin_pred)
lin_error = mean_squared_error(y_test, lin_pred)

print('Accuracy (Linear): ', lin_accuracy)
print('F1 (Linear): ', lin_f1)
print('Mean squared error: ', lin_error)

Accuracy (Polynomial Kernel):  0.475
F1 (Polynomial Kernel):  0.6440677966101694
Mean squared error:  0.525
---------
Accuracy (SVC standard):  0.975
F1 (SVC standard):  0.972972972972973
Mean squared error standard:  0.025
---------
Accuracy (Linear):  0.975
F1 (Linear):  0.9743589743589743
Mean squared error:  0.025


df = pd.read_csv(r"C:\Users\marti\OneDrive\AI_Lab\HCC1806_SmartS_Filtered_Normalised_3000_Data_train.txt", delimiter="\ ",engine='python',index_col=0)
df = name_clean(df)
df = df.T


H_Smart_meta = pd.read_csv(r"C:\Users\marti\OneDrive\AI_Lab\HCC1806_SmartS_MetaData.tsv", delimiter="\t",engine='python',index_col=0)
H_Smart_meta = name_clean(H_Smart_meta.T)
H_Smart_meta = H_Smart_meta.drop(index= ["PreprocessingTag", "ProcessingComments", "Hours", "PCR Plate", "Cell Line","Pos"])


# Create a dictionary to map cell names to their condition
condition_dict = H_Smart_meta.T['Condition'].to_dict()

# Use the `map` function to add the condition as a new column in df
df['Condition'] = df.index.map(condition_dict)
df['Condition'] = df['Condition'].map({'Hypo': 1, 'Normo': 0})


X = df.drop(columns=['Condition'])
y = df['Condition']


X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.2, random_state=42)


colors = ['#04395E', '#C3C49E']
fig, ax = plt.subplots(figsize=(8, 4))
sns.countplot(x = 'Condition', data = df, palette = 'viridis')
plt.show()
df['Condition'].value_counts()

1    97
0    85
Name: Condition, dtype: int64


wcss = []
for i in range(1, 10):
    kmeans = KMeans(n_clusters=i, init='k-means++', max_iter=300, n_init=10, random_state=0)
    kmeans.fit(df.drop(columns=['Condition']))
    wcss.append(kmeans.inertia_)
plt.plot(range(1, 10), wcss, color = 'purple')
plt.title('Elbow Method')
plt.axvline(x = 2, color='forestgreen', linestyle='--', alpha = 0.7)
plt.xlabel('Number of clusters')
plt.ylabel('WCSS')
plt.show()


kmeans = KMeans(2, init='k-means++',  n_init=2000, random_state=42).fit(X)


cm = confusion_matrix(df['Condition'], kmeans.labels_)
disp = ConfusionMatrixDisplay(confusion_matrix=cm)
disp.plot()

<sklearn.metrics._plot.confusion_matrix.ConfusionMatrixDisplay at 0x18cce184c70>


accuracy_score(kmeans.labels_,y)

0.4835164835164835


hierarchical = AgglomerativeClustering(
                                distance_threshold = 0.15, 
                                n_clusters=None)
hierarchical.fit(X)
Z = linkage(X, 'ward')

plt.figure(figsize=(10, 7), facecolor='white')
plt.style.use('default')
plt.title("Dendrogram")
dendrogram(Z)
plt.gca().set_facecolor('white')

plt.show()


non_tuned_model_rf = RandomForestClassifier()
non_tuned_model_rf.fit(X_train, y_train)

RandomForestClassifier()

RandomForestClassifier()


important_features = pd.DataFrame(non_tuned_model_rf.feature_importances_, index = df.drop(columns=['Condition']).columns).rename(columns = {0: 'Important features'})
important_features = important_features[important_features['Important features'] > 0]
important_features


featured_df = df[important_features.T.columns]


# Use the `map` function to add the condition as a new column in df
featured_df['Condition'] = df.index.map(condition_dict)
featured_df['Condition'] = featured_df['Condition'].map({'Hypo': 1, 'Normo': 0})


X = featured_df.drop(columns=['Condition'])
y = featured_df['Condition']
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.2, random_state=29)


non_tuned_model_lr = LogisticRegression(random_state=29)
non_tuned_model_lr.fit(X_train, y_train)
print("Non-tuned logistic score score on train data:", non_tuned_model_lr.score(X_train,y_train))
print("Non-tuned logistic score on validation data:",non_tuned_model_lr.score(X_test,y_test))
print("---")
non_tuned_model_rf = RandomForestClassifier(random_state=29)
non_tuned_model_rf.fit(X_train, y_train)
print("Non-tuned random forest score on train data:", non_tuned_model_rf.score(X_train,y_train))
print("Non-tuned random forest score on validation data:",non_tuned_model_rf.score(X_test,y_test))
print("---")
non_tuned_model_svc = svm.LinearSVC(random_state=29)
non_tuned_model_svc.fit(X_train,y_train)
print("Non-tuned Linear SVC score on train data:", non_tuned_model_svc.score(X_train,y_train))
print("Non-tuned Linear SVC score on validation data:",non_tuned_model_svc.score(X_test,y_test))

Non-tuned logistic score score on train data: 1.0
Non-tuned logistic score on validation data: 0.972972972972973
---
Non-tuned random forest score on train data: 1.0
Non-tuned random forest score on validation data: 0.972972972972973
---
Non-tuned Linear SVC score on train data: 1.0
Non-tuned Linear SVC score on validation data: 0.972972972972973


model = Sequential([
    Dense(64, activation='relu', input_dim=X_train.shape[1]),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(1, activation='sigmoid'),
])
model.compile(optimizer='Nadam', loss='binary_crossentropy', metrics=['binary_accuracy'])
model.fit(X_train, y_train, epochs=30, batch_size=1, validation_data=(X_test, y_test))
loss, accuracy = model.evaluate(X_test, y_test)
print(f'Neural network score on validation data: {accuracy}')

Epoch 1/30
145/145 [==============================] - 3s 3ms/step - loss: 34.9595 - binary_accuracy: 0.8552 - val_loss: 1.4756 - val_binary_accuracy: 0.9730
Epoch 2/30
145/145 [==============================] - 0s 2ms/step - loss: 8.3203 - binary_accuracy: 0.9448 - val_loss: 0.0032 - val_binary_accuracy: 1.0000
Epoch 3/30
145/145 [==============================] - 0s 1ms/step - loss: 1.9659 - binary_accuracy: 0.9793 - val_loss: 17.7709 - val_binary_accuracy: 0.9730
Epoch 4/30
145/145 [==============================] - 0s 1ms/step - loss: 1.4260e-14 - binary_accuracy: 1.0000 - val_loss: 16.9214 - val_binary_accuracy: 0.9730
Epoch 5/30
145/145 [==============================] - 0s 1ms/step - loss: 2.0413e-14 - binary_accuracy: 1.0000 - val_loss: 16.9235 - val_binary_accuracy: 0.9730
Epoch 6/30
145/145 [==============================] - 0s 1ms/step - loss: 1.7672e-14 - binary_accuracy: 1.0000 - val_loss: 16.9253 - val_binary_accuracy: 0.9730
Epoch 7/30
145/145 [==============================] - 0s 1ms/step - loss: 1.5607e-14 - binary_accuracy: 1.0000 - val_loss: 16.9270 - val_binary_accuracy: 0.9730
Epoch 8/30
145/145 [==============================] - 0s 1ms/step - loss: 1.3978e-14 - binary_accuracy: 1.0000 - val_loss: 16.9284 - val_binary_accuracy: 0.9730
Epoch 9/30
145/145 [==============================] - 0s 1ms/step - loss: 1.2668e-14 - binary_accuracy: 1.0000 - val_loss: 16.9296 - val_binary_accuracy: 0.9730
Epoch 10/30
145/145 [==============================] - 0s 1ms/step - loss: 1.1665e-14 - binary_accuracy: 1.0000 - val_loss: 16.9310 - val_binary_accuracy: 0.9730
Epoch 11/30
145/145 [==============================] - 0s 1ms/step - loss: 1.0669e-14 - binary_accuracy: 1.0000 - val_loss: 16.9319 - val_binary_accuracy: 0.9730
Epoch 12/30
145/145 [==============================] - 0s 1ms/step - loss: 9.8987e-15 - binary_accuracy: 1.0000 - val_loss: 16.9330 - val_binary_accuracy: 0.9730
Epoch 13/30
145/145 [==============================] - 0s 1ms/step - loss: 9.2316e-15 - binary_accuracy: 1.0000 - val_loss: 16.9341 - val_binary_accuracy: 0.9730
Epoch 14/30
145/145 [==============================] - 0s 1ms/step - loss: 8.6517e-15 - binary_accuracy: 1.0000 - val_loss: 16.9349 - val_binary_accuracy: 0.9730
Epoch 15/30
145/145 [==============================] - 0s 1ms/step - loss: 8.1437e-15 - binary_accuracy: 1.0000 - val_loss: 16.9358 - val_binary_accuracy: 0.9730
Epoch 16/30
145/145 [==============================] - 0s 1ms/step - loss: 7.6878e-15 - binary_accuracy: 1.0000 - val_loss: 16.9366 - val_binary_accuracy: 0.9730
Epoch 17/30
145/145 [==============================] - 0s 2ms/step - loss: 7.2842e-15 - binary_accuracy: 1.0000 - val_loss: 16.9373 - val_binary_accuracy: 0.9730
Epoch 18/30
145/145 [==============================] - 0s 1ms/step - loss: 6.9226e-15 - binary_accuracy: 1.0000 - val_loss: 16.9381 - val_binary_accuracy: 0.9730
Epoch 19/30
145/145 [==============================] - 0s 1ms/step - loss: 6.5933e-15 - binary_accuracy: 1.0000 - val_loss: 16.9388 - val_binary_accuracy: 0.9730
Epoch 20/30
145/145 [==============================] - 0s 1ms/step - loss: 6.2963e-15 - binary_accuracy: 1.0000 - val_loss: 16.9394 - val_binary_accuracy: 0.9730
Epoch 21/30
145/145 [==============================] - 0s 1ms/step - loss: 6.0247e-15 - binary_accuracy: 1.0000 - val_loss: 16.9400 - val_binary_accuracy: 0.9730
Epoch 22/30
145/145 [==============================] - 0s 1ms/step - loss: 5.7759e-15 - binary_accuracy: 1.0000 - val_loss: 16.9405 - val_binary_accuracy: 0.9730
Epoch 23/30
145/145 [==============================] - 0s 1ms/step - loss: 5.5481e-15 - binary_accuracy: 1.0000 - val_loss: 16.9412 - val_binary_accuracy: 0.9730
Epoch 24/30
145/145 [==============================] - 0s 1ms/step - loss: 5.3367e-15 - binary_accuracy: 1.0000 - val_loss: 16.9415 - val_binary_accuracy: 0.9730
Epoch 25/30
145/145 [==============================] - 0s 1ms/step - loss: 5.1414e-15 - binary_accuracy: 1.0000 - val_loss: 16.9420 - val_binary_accuracy: 0.9730
Epoch 26/30
145/145 [==============================] - 0s 1ms/step - loss: 4.9595e-15 - binary_accuracy: 1.0000 - val_loss: 16.9426 - val_binary_accuracy: 0.9730
Epoch 27/30
145/145 [==============================] - 0s 1ms/step - loss: 4.7911e-15 - binary_accuracy: 1.0000 - val_loss: 16.9430 - val_binary_accuracy: 0.9730
Epoch 28/30
145/145 [==============================] - 0s 1ms/step - loss: 4.6865e-15 - binary_accuracy: 1.0000 - val_loss: 16.9438 - val_binary_accuracy: 0.9730
Epoch 29/30
145/145 [==============================] - 0s 1ms/step - loss: 4.4817e-15 - binary_accuracy: 1.0000 - val_loss: 16.9442 - val_binary_accuracy: 0.9730
Epoch 30/30
145/145 [==============================] - 0s 1ms/step - loss: 4.3448e-15 - binary_accuracy: 1.0000 - val_loss: 16.9447 - val_binary_accuracy: 0.9730
2/2 [==============================] - 0s 4ms/step - loss: 16.9447 - binary_accuracy: 0.9730
Neural network score on validation data: 0.9729729890823364


# Accuracy scores
lr_train_accuracy = accuracy_score(y_train, non_tuned_model_lr.predict(X_train))
lr_val_accuracy = accuracy_score(y_test, non_tuned_model_lr.predict(X_test))
rf_train_accuracy = accuracy_score(y_train, non_tuned_model_rf.predict(X_train))
rf_val_accuracy = accuracy_score(y_test, non_tuned_model_rf.predict(X_test))
svc_train_accuracy = accuracy_score(y_train, non_tuned_model_svc.predict(X_train))
svc_val_accuracy = accuracy_score(y_test, non_tuned_model_svc.predict(X_test))
nn_train_accuracy = accuracy_score(y_train, np.round(model.predict(X_train, verbose=0)))
nn_val_accuracy = accuracy_score(y_test, np.round(model.predict(X_test, verbose=0)))

# Precision scores
lr_train_precision = precision_score(y_train, non_tuned_model_lr.predict(X_train))
lr_val_precision = precision_score(y_test, non_tuned_model_lr.predict(X_test))
rf_train_precision = precision_score(y_train, non_tuned_model_rf.predict(X_train))
rf_val_precision = precision_score(y_test, non_tuned_model_rf.predict(X_test))
svc_train_precision = precision_score(y_train, non_tuned_model_svc.predict(X_train))
svc_val_precision = precision_score(y_test, non_tuned_model_svc.predict(X_test))
nn_train_precision = precision_score(y_train, np.round(model.predict(X_train)))
nn_val_precision = precision_score(y_test, np.round(model.predict(X_test)))

# Recall scores
lr_train_recall = recall_score(y_train, non_tuned_model_lr.predict(X_train))
lr_val_recall = recall_score(y_test, non_tuned_model_lr.predict(X_test))
rf_train_recall = recall_score(y_train, non_tuned_model_rf.predict(X_train))
rf_val_recall = recall_score(y_test, non_tuned_model_rf.predict(X_test))
svc_train_recall = recall_score(y_train, non_tuned_model_svc.predict(X_train))
svc_val_recall = recall_score(y_test, non_tuned_model_svc.predict(X_test))
nn_train_recall = recall_score(y_train, np.round(model.predict(X_train, verbose=0)))
nn_val_recall = recall_score(y_test, np.round(model.predict(X_test, verbose=0)))

# F1 scores
lr_train_f1 = f1_score(y_train, non_tuned_model_lr.predict(X_train))
lr_val_f1 = f1_score(y_test, non_tuned_model_lr.predict(X_test))
rf_train_f1 = f1_score(y_train, non_tuned_model_rf.predict(X_train))
rf_val_f1 = f1_score(y_test, non_tuned_model_rf.predict(X_test))
svc_train_f1 = f1_score(y_train, non_tuned_model_svc.predict(X_train))
svc_val_f1 = f1_score(y_test, non_tuned_model_svc.predict(X_test))
nn_train_f1 = f1_score(y_train, np.round(model.predict(X_train, verbose=0)))
nn_val_f1 = f1_score(y_test, np.round(model.predict(X_test, verbose=0)))

# AUC scores
lr_train_auc = metrics.roc_auc_score(y_train, non_tuned_model_lr.predict(X_train))
lr_val_auc = metrics.roc_auc_score(y_test, non_tuned_model_lr.predict(X_test))
rf_train_auc = metrics.roc_auc_score(y_train, non_tuned_model_rf.predict(X_train))
rf_val_auc = metrics.roc_auc_score(y_test, non_tuned_model_rf.predict(X_test))
svc_train_auc = metrics.roc_auc_score(y_train, non_tuned_model_svc.predict(X_train))
svc_val_auc = metrics.roc_auc_score(y_test, non_tuned_model_svc.predict(X_test))
nn_train_auc = metrics.roc_auc_score(y_train, np.round(model.predict(X_train, verbose=0)))
nn_val_auc = metrics.roc_auc_score(y_test, np.round(model.predict(X_test, verbose=0)))

# Define the labels   
score_labels = ['LR on Train', 'LR on Validation', 
                'RF on Train', 'RF on Validation', 
                'SVC on Train', 'SVC on Validation',
                'NN on Train', 'NN on Validation']

# Define scores
accuracy_scores = [lr_train_accuracy, lr_val_accuracy,
                   rf_train_accuracy, rf_val_accuracy,
                   svc_train_accuracy, svc_val_accuracy,
                   nn_train_accuracy, nn_val_accuracy]

precision_scores = [lr_train_precision, lr_val_precision,
                    rf_train_precision, rf_val_precision,
                    svc_train_precision, svc_val_precision,
                    nn_train_precision, nn_val_precision]

recall_scores = [lr_train_recall, lr_val_recall,
                 rf_train_recall, rf_val_recall,
                 svc_train_recall, svc_val_recall,
                 nn_train_recall, nn_val_recall]

f1_scores = [lr_train_f1, lr_val_f1,
             rf_train_f1, rf_val_f1,
             svc_train_f1, svc_val_f1,
             nn_train_f1, nn_val_f1]

auc_scores = [lr_train_auc, lr_val_auc,
              rf_train_auc, rf_val_auc,
              svc_train_auc, svc_val_auc,
              nn_train_auc, nn_val_auc]

# Create a grouped bar plot
bar_width = 0.12
index = np.arange(len(score_labels))  # Adjust the index to account for number of models

plt.figure(figsize=(15, 7))
plt.bar(index, accuracy_scores, width=bar_width, label='Accuracy', alpha = 0.8, color = 'indigo')
plt.bar(index + bar_width, precision_scores, width=bar_width, label='Precision', alpha = 0.8, color = 'darkslateblue')
plt.bar(index + 2 * bar_width, recall_scores, width=bar_width, label='Recall',alpha = 0.8, color = 'teal')
plt.bar(index + 3 * bar_width, f1_scores, width=bar_width, label='F1 Score', alpha = 0.8,color = 'seagreen')
plt.bar(index + 4 * bar_width, auc_scores, width=bar_width, label='AUC', alpha = 0.8, color = 'gold')

plt.xlabel('Models')
plt.ylabel('Scores')
plt.title('Model Evaluation Scores')
plt.xticks(index + 2 * bar_width, score_labels, rotation = 45)  # Rotate xticks for better readability
plt.legend()
plt.tight_layout()
plt.gca().set_axisbelow(True)
plt.grid(True)
plt.show()

5/5 [==============================] - 0s 983us/step
2/2 [==============================] - 0s 2ms/step


print("Metrics on Logistic Regression")
print("Accuracy score on validation data: ", lr_val_accuracy)
print("Precision score on validation data: ", lr_val_precision)
print("Recall score on validation data: ", lr_val_recall)
print("F1 score on validation data: ", lr_val_f1)
print("AUC on validation data: ", lr_val_auc)
cross_val = cross_val_score(non_tuned_model_lr, X_train, y_train, scoring = 'accuracy', cv = 20)
print("Minimum  cross validation score:", min(cross_val)) 
print("Maximum cross validation score:", max(cross_val))
print("Average cross validation score:", np.mean(cross_val))

Metrics on Logistic Regression
Accuracy score on validation data:  0.972972972972973
Precision score on validation data:  1.0
Recall score on validation data:  0.9444444444444444
F1 score on validation data:  0.9714285714285714
AUC on validation data:  0.9722222222222222
Minimum  cross validation score: 1.0
Maximum cross validation score: 1.0
Average cross validation score: 1.0


M_Smart_norm = pd.read_csv(r"C:\Users\marti\OneDrive\AI_Lab\MCF7_SmartS_Filtered_Normalised_3000_Data_train.txt",
                           delimiter = "\ ", engine="python", index_col=0)
M_Smart_metadata =pd.read_csv(r"C:\Users\marti\OneDrive\AI_Lab\MCF7_SmartS_MetaData.tsv",
                              delimiter = "\t", engine="python", index_col=0)


M_Smart_norm = name_clean(M_Smart_norm)
M_Smart_norm.head()


M_Smart_metadata = name_clean(M_Smart_metadata.T).T
M_Smart_metadata.head()


label = pd.Series(np.zeros(383), index = M_Smart_metadata.index)
for i in M_Smart_metadata.index:
    if M_Smart_metadata.Condition[i] == 'Norm':
        label[i] = 0
    elif M_Smart_metadata.Condition[i] == 'Hypo':
        label[i] = 1
label = label.astype('int64')
label

Filename
S28     1
S29     1
S30     1
S1      0
S2      0
       ..
S359    0
S360    0
S379    1
S380    1
S381    1
Length: 383, dtype: int64


M_Smart_norm_complete = M_Smart_norm.T.copy()
M_Smart_norm_complete['Condition'] = label
M_Smart_norm_complete = M_Smart_norm_complete.T


pca2 = PCA(n_components=2)
M_Smart_norm_2 = pca2.fit_transform(M_Smart_norm.T)
colors = ['purple', 'forestgreen']
labels = [0, 1]
plt.scatter(M_Smart_norm_2[:, 0], M_Smart_norm_2[:, 1], 
            c=[colors[label] for label in M_Smart_norm_complete.T.Condition])
legend_handles = [mpatches.Patch(color=color, label=label) for color, label in zip(colors, labels)]
plt.legend(handles=legend_handles)
plt.title('Scatter plot')
plt.xlabel('PCA1')
plt.ylabel('PCA2')
plt.show()


Xtrain, Xtest, ytrain, ytest = model_selection.train_test_split(M_Smart_norm.T, M_Smart_norm_complete.T.Condition, test_size=0.3, random_state=42)


clf = RandomForestClassifier(random_state=43)
clf.fit(Xtrain,ytrain)

RandomForestClassifier(random_state=43)

RandomForestClassifier(random_state=43)


important_features = pd.DataFrame(clf.feature_importances_, index = M_Smart_norm.T.columns).rename(columns = {0: 'Important features'})
important_features = important_features[important_features['Important features'] > 0]


important_features


M_Smart_norm_important = M_Smart_norm.T.copy()
M_Smart_norm_important = M_Smart_norm_important.loc[:, important_features.index.tolist()]
M_Smart_norm_important_complete = M_Smart_norm_important.copy()
M_Smart_norm_important_complete['Condition'] = label


Xtrain, Xtest, ytrain, ytest = model_selection.train_test_split(M_Smart_norm_important, 
                                                M_Smart_norm_important_complete.Condition, test_size=0.3, random_state=42)


param_grid_1= [{'penalty': ['l2', None], 'C':list(np.logspace(-4, 4, 10)) + [1] ,
                'solver' : ['lbfgs', 'newton-cg', 'sag', 'saga'], 
                       'max_iter': [200,800,1500,4000], 'fit_intercept' : [True, False] }]
param_grid_2 = [{'penalty': ['l2','l1'], 'C':np.logspace(-4, 4, 10) + [1] , 'solver' : ['liblinear'], 
                       'max_iter': [200,800,1500,4000], 'fit_intercept' : [True, False]}]


#Empty model
log_reg = LogisticRegression(random_state = 42)
log_reg.fit(Xtrain, ytrain)
# First Parameter Grid
grid_1 = RandomizedSearchCV(log_reg, param_grid_1, cv = 5, scoring = 'accuracy', n_jobs = -1)
grid_1.fit(Xtrain, ytrain)
# Second Parameter Grid
grid_2 = RandomizedSearchCV(log_reg, param_grid_2, cv = 5, scoring = 'accuracy', n_jobs = -1)
grid_2.fit(Xtrain, ytrain)

RandomizedSearchCV(cv=5, estimator=LogisticRegression(random_state=42),
                   n_jobs=-1,
                   param_distributions=[{'C': array([1.00010000e+00, 1.00077426e+00, 1.00599484e+00, 1.04641589e+00,
       1.35938137e+00, 3.78255940e+00, 2.25443469e+01, 1.67810054e+02,
       1.29254967e+03, 1.00010000e+04]),
                                         'fit_intercept': [True, False],
                                         'max_iter': [200, 800, 1500, 4000],
                                         'penalty': ['l2', 'l1'],
                                         'solver': ['liblinear']}],
                   scoring='accuracy')

RandomizedSearchCV(cv=5, estimator=LogisticRegression(random_state=42),
                   n_jobs=-1,
                   param_distributions=[{'C': array([1.00010000e+00, 1.00077426e+00, 1.00599484e+00, 1.04641589e+00,
       1.35938137e+00, 3.78255940e+00, 2.25443469e+01, 1.67810054e+02,
       1.29254967e+03, 1.00010000e+04]),
                                         'fit_intercept': [True, False],
                                         'max_iter': [200, 800, 1500, 4000],
                                         'penalty': ['l2', 'l1'],
                                         'solver': ['liblinear']}],
                   scoring='accuracy')

LogisticRegression(random_state=42)

LogisticRegression(random_state=42)


fig, axes = plt.subplots(1, 3, figsize=(12, 5))

fig.suptitle('Confusion Matrix on Train data', fontweight = 'bold')
# Confusion Matrix and score of Empty Model
y_train_pred_1 = model_selection.cross_val_predict(log_reg, Xtrain, ytrain, cv = 5)
empty_model_cm = ConfusionMatrixDisplay.from_predictions(ytrain, y_train_pred_1, ax = axes[0], colorbar=False)
axes[0].set_title('Empty Model')
axes[0].set_xlabel('Predicted label')
axes[0].set_ylabel('True label')

# Confusion Matrix and score of First Grid model
y_train_pred_2 = model_selection.cross_val_predict(grid_1, Xtrain, ytrain, cv = 5)
first_grid_cm =ConfusionMatrixDisplay.from_predictions(ytrain, y_train_pred_2, ax = axes[1],  colorbar=False)
axes[1].set_title('First Grid Model')
axes[1].set_xlabel('Predicted label')
axes[1].set_ylabel('True label')

# Confusion Matrix and score of Second Grid model
y_train_pred_3 = model_selection.cross_val_predict(grid_2, Xtrain, ytrain, cv = 5)
second_grid_cm = ConfusionMatrixDisplay.from_predictions(ytrain, y_train_pred_3, ax = axes[2],  colorbar=False)
axes[2].set_title('Second Grid Model')
axes[2].set_xlabel('Predicted label')
axes[2].set_ylabel('True label')

plt.tight_layout()
plt.show()


# Test Performance on empty model
y_test_pred_1 = log_reg.predict(Xtest)
print("\033[1m"+"Empty model"+ "\033[0m")
print('The accuracy score in empty model is:', metrics.accuracy_score(ytest, y_test_pred_1))
print('The precision score in empty model is:', metrics.precision_score(ytest, y_test_pred_1))
print('The recall score in empty model is: ',metrics.recall_score(ytest, y_test_pred_1))
print('The f1 score in empty model is: ',metrics.f1_score(ytest, y_test_pred_1))
print('The roc auc score in empty model is', metrics.roc_auc_score(ytest, y_test_pred_1))

# Test Performance on First Grid model
y_test_pred_2 = grid_1.predict(Xtest)
print("\n\033[1m"+"First Grid model" + "\033[0m")
print('The accuracy score in first grid is:', metrics.accuracy_score(ytest, y_test_pred_2))
print('The precision score in first grid is:', metrics.precision_score(ytest, y_test_pred_2))
print('The recall score in first grid is: ',metrics.recall_score(ytest, y_test_pred_2))
print('The f1 score in first grid is: ',metrics.f1_score(ytest, y_test_pred_2))
print('The roc auc score in first grid is', metrics.roc_auc_score(ytest, y_test_pred_2))

# Test Performance on Second Grid model
y_test_pred_3 = grid_2.predict(Xtest)
print("\n\033[1m"+"Second Grid model" + "\033[0m")
print('The accuracy score in second grid is:', metrics.accuracy_score(ytest, y_test_pred_3))
print('The precision score in second grid is:', metrics.precision_score(ytest, y_test_pred_3))
print('The recall score in second grid is: ',metrics.recall_score(ytest, y_test_pred_3))
print('The f1 score in second grid is: ',metrics.f1_score(ytest, y_test_pred_3))
print('The roc auc score in second grid is', metrics.roc_auc_score(ytest, y_test_pred_3))

Empty model
The accuracy score in empty model is: 1.0
The precision score in empty model is: 1.0
The recall score in empty model is:  1.0
The f1 score in empty model is:  1.0
The roc auc score in empty model is 1.0

First Grid model
The accuracy score in first grid is: 1.0
The precision score in first grid is: 1.0
The recall score in first grid is:  1.0
The f1 score in first grid is:  1.0
The roc auc score in first grid is 1.0

Second Grid model
The accuracy score in second grid is: 1.0
The precision score in second grid is: 1.0
The recall score in second grid is:  1.0
The f1 score in second grid is:  1.0
The roc auc score in second grid is 1.0


H_drop_fil = pd.read_csv(r'C:\Users\marti\OneDrive\AI_Lab\HCC1806_DropS_Filtered_Normalised_3000_Data_train.txt', sep = ' ')

H_drop_fil.head()


print('Dimension of the dataset: ', H_drop_fil.shape)
n_rows= H_drop_fil.shape[0]
n_col= H_drop_fil.shape[1]
print('Number of rows:', n_rows, '\nNumber of columns:', n_col)

Dimension of the dataset:  (3000, 14682)
Number of rows: 3000 
Number of columns: 14682


if len(H_drop_fil.dtypes[H_drop_fil.dtypes == 'int64']) == n_col:
    print('All entries are integer numbers')

def empty_cols(df):
    nulldf=[]
    for col in df.columns:
        if df[col].sum()==0:
            nulldf.append(col)
    return nulldf
print('Empty columns:', empty_cols(H_drop_fil))

All entries are integer numbers
Empty columns: []


print(H_drop_fil.isna().sum().sum())
print(H_drop_fil.isnull().sum().sum())

0
0


duplicate_cols = []
for i in H_drop_fil.duplicated():
    if i == True:
        duplicate_cols.append(i)
print('Duplicate columns: ' +str(duplicate_cols))
duplicate_rows = []
for i in H_drop_fil.T.duplicated():
    if i == True:
        duplicate_rows.append(i)
print(f"Duplicate rows: {duplicate_rows}")

Duplicate columns: []
Duplicate rows: []


H_drop_fil.describe()


min_mean= min(H_drop_fil.describe().loc['mean'].tolist())
max_mean= max(H_drop_fil.describe().loc['mean'].tolist())
print(f"The mean ranges from {min_mean} to {max_mean}")
std= np.mean(np.array(H_drop_fil.describe().loc['std'].tolist()))
print(f"The average standard deviation is {std}")

extr_low= min(H_drop_fil.describe().loc['min'].tolist())
extr_high= max(H_drop_fil.describe().loc['max'].tolist())
print(f"The extreme values of the dataset are {extr_low} and {extr_high}")

The mean ranges from 0.013 to 0.07333333333333333
The average standard deviation is 0.2790719033977203
The extreme values of the dataset are 0.0 and 110.0


plt.figure(figsize=(17, 12))
corr = H_drop_fil[0:50].corr()
sns.heatmap(corr, cmap = 'viridis');
plt.title('Correlation matrix');


df = pd.DataFrame(H_drop_fil.T)
df['label'] = np.zeros(H_drop_fil.shape[1])
for i in H_drop_fil:
    if 'Normoxia' in i:
        df['label'][i] = 0
    if 'Hypoxia' in i:
        df['label'][i] = 1

df['label'] = df['label'].astype(int)
df = df.T


plt.figure()
sns.countplot(data= df.T, x ='label',  palette = ['forestgreen', 'purple'], hue = 'label');
plt.title("Distribution of Normoxic and Hypoxic Cells");
plt.legend(title='Condition', loc='best', labels=['Normo', 'Hypo']);
plt.xticks([]);


np.bincount(df.T['label'])

array([5783, 8899], dtype=int64)


normo = df.T[df.T['label'] == 0].T
normo.drop('label', axis = 'index', inplace = True)
hypo = df.T[df.T['label'] == 1].T
hypo.drop('label', axis = 'index', inplace = True)


sns.kdeplot(data=normo.mean(), label = 'Normo', clip = [0,1000], color = 'purple')
sns.kdeplot(data=hypo.mean(), label = 'Hypo', clip = [0,1000], color = 'forestgreen')
plt.legend()
plt.show()


# Perform PCA on the unfiltered dataset
pca = PCA(n_components=2)
pca_result = pca.fit_transform(df.T)

# Get condition labels for the samples
conditions = df.T["label"]

# Create a PCA plot
plt.figure(figsize=(8, 6))
for condition, color in zip(conditions.unique(), ['green', 'purple']):
    idx = conditions == condition
    plt.scatter(pca_result[idx, 0], pca_result[idx, 1], c=color, label=condition, alpha = 0.3)
plt.legend()
plt.xlabel('PC1')
plt.ylabel('PC2')
plt.title('2D plot')
plt.show()


# PCA model with three components
import plotly.express as px

condition = df.T['label'].astype(str)
pca = PCA(n_components=3)
pca_result = pca.fit_transform(df.T)
g = pd.DataFrame(pca_result)
g.set_index(condition, inplace =True)


fig = px.scatter_3d(g, x=0, y=1, z=2,
                     color = g.index, color_discrete_sequence = ['purple', 'green']
                    
                    )
#
fig.show()


X = df.T.drop('label', axis = 1)
y = df.T['label']
X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size = 0.25, random_state=43)


clf = RandomForestClassifier(random_state=43)
clf.fit(X_train,y_train)

RandomForestClassifier(random_state=43)

RandomForestClassifier(random_state=43)


important_features = pd.DataFrame(clf.feature_importances_, index = df.T.drop('label', axis = 1).columns).rename(columns = {0: 'Important features'})
important_features = important_features[important_features['Important features'] > 0]
important_features


X = X[important_features.index]


model = Sequential([
    Dense(64, activation='relu', input_dim=X_train.shape[1]),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(1, activation='sigmoid'),
])
model.compile(optimizer='Nadam', loss='binary_crossentropy', metrics=['binary_accuracy'])
model.fit(X_train, y_train, epochs=10, batch_size=1, validation_data=(X_test, y_test))
loss, accuracy = model.evaluate(X_test, y_test)
print(f'Test accuracy: {accuracy}')

Epoch 1/10
11011/11011 [==============================] - 103s 9ms/step - loss: 0.1721 - binary_accuracy: 0.9355 - val_loss: 0.1268 - val_binary_accuracy: 0.9504
Epoch 2/10
11011/11011 [==============================] - 98s 9ms/step - loss: 0.0874 - binary_accuracy: 0.9697 - val_loss: 0.1163 - val_binary_accuracy: 0.9600
Epoch 3/10
11011/11011 [==============================] - 122s 11ms/step - loss: 0.0468 - binary_accuracy: 0.9845 - val_loss: 0.1636 - val_binary_accuracy: 0.9567
Epoch 4/10
11011/11011 [==============================] - 125s 11ms/step - loss: 0.0249 - binary_accuracy: 0.9918 - val_loss: 0.2268 - val_binary_accuracy: 0.9589
Epoch 5/10
11011/11011 [==============================] - 119s 11ms/step - loss: 0.0173 - binary_accuracy: 0.9965 - val_loss: 0.4105 - val_binary_accuracy: 0.9548
Epoch 6/10
11011/11011 [==============================] - 123s 11ms/step - loss: 0.0135 - binary_accuracy: 0.9968 - val_loss: 0.4535 - val_binary_accuracy: 0.9594
Epoch 7/10
11011/11011 [==============================] - 122s 11ms/step - loss: 0.0081 - binary_accuracy: 0.9982 - val_loss: 0.8072 - val_binary_accuracy: 0.9521
Epoch 8/10
11011/11011 [==============================] - 121s 11ms/step - loss: 0.0086 - binary_accuracy: 0.9976 - val_loss: 1.0051 - val_binary_accuracy: 0.9534
Epoch 9/10
11011/11011 [==============================] - 122s 11ms/step - loss: 0.0072 - binary_accuracy: 0.9987 - val_loss: 0.8627 - val_binary_accuracy: 0.9523
Epoch 10/10
11011/11011 [==============================] - 124s 11ms/step - loss: 0.0121 - binary_accuracy: 0.9981 - val_loss: 0.5415 - val_binary_accuracy: 0.9567
115/115 [==============================] - 1s 5ms/step - loss: 0.5415 - binary_accuracy: 0.9567
Test accuracy: 0.9566875696182251


y_pred1 = model.predict(X_test)
print('Accuracy Score: ' + str(accuracy))
print('Binary Cross Entropy  loss: ' +str(loss))
print('Mean Squared Error Loss: ' + str(mean_squared_error(y_test, y_pred1)))

115/115 [==============================] - 0s 3ms/step
Accuracy Score: 0.9566875696182251
Binary Cross Entropy  loss: 0.541532039642334
Mean Squared Error Loss: 0.03805433755640302


def create_model(optimizer):
    # define model
    model = Sequential([
    Dense(64, activation='relu', input_dim=X_train.shape[1]),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(1, activation='sigmoid'),
])
    model.compile(optimizer='Nadam', loss='binary_crossentropy', metrics=['binary_accuracy'])
    return model


model = KerasClassifier(model=create_model, loss="binary_crossentropy", epochs=10, batch_size=10, verbose=0)
optimizer = ['SGD', 'RMSprop', 'Adagrad', 'Adadelta', 'Adam', 'Adamax', 'Nadam']
param_grid = dict(model__optimizer=optimizer)
grid = GridSearchCV(estimator=model, param_grid=param_grid, n_jobs=-1, cv=5)
grid_result = grid.fit(X_train, y_train)
print(grid_result.best_params_)

{'model__optimizer': 'Adagrad'}


X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size=0.25, random_state=43)


model = Sequential([
    Dense(64, activation='relu', input_dim=X_train.shape[1]),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(1, activation='sigmoid'),
])
model.compile(optimizer='Adagrad', loss='binary_crossentropy', metrics=['binary_accuracy'])
model.fit(X_train, y_train, epochs=30, batch_size=1, validation_data=(X_test, y_test))
loss, accuracy = model.evaluate(X_test, y_test)
print(f'Test accuracy: {accuracy}')

Epoch 1/30
11011/11011 [==============================] - 10s 845us/step - loss: 0.2820 - binary_accuracy: 0.8859 - val_loss: 0.1563 - val_binary_accuracy: 0.9412
Epoch 2/30
11011/11011 [==============================] - 9s 815us/step - loss: 0.1258 - binary_accuracy: 0.9531 - val_loss: 0.1305 - val_binary_accuracy: 0.9480
Epoch 3/30
11011/11011 [==============================] - 9s 859us/step - loss: 0.0992 - binary_accuracy: 0.9647 - val_loss: 0.1179 - val_binary_accuracy: 0.9515
Epoch 4/30
11011/11011 [==============================] - 11s 1ms/step - loss: 0.0841 - binary_accuracy: 0.9704 - val_loss: 0.1113 - val_binary_accuracy: 0.9559
Epoch 5/30
11011/11011 [==============================] - 9s 815us/step - loss: 0.0734 - binary_accuracy: 0.9745 - val_loss: 0.1079 - val_binary_accuracy: 0.9583
Epoch 6/30
11011/11011 [==============================] - 9s 842us/step - loss: 0.0654 - binary_accuracy: 0.9780 - val_loss: 0.1056 - val_binary_accuracy: 0.9597
Epoch 7/30
11011/11011 [==============================] - 10s 863us/step - loss: 0.0584 - binary_accuracy: 0.9817 - val_loss: 0.1046 - val_binary_accuracy: 0.9602
Epoch 8/30
11011/11011 [==============================] - 9s 822us/step - loss: 0.0531 - binary_accuracy: 0.9837 - val_loss: 0.1036 - val_binary_accuracy: 0.9619
Epoch 9/30
11011/11011 [==============================] - 9s 822us/step - loss: 0.0481 - binary_accuracy: 0.9863 - val_loss: 0.1034 - val_binary_accuracy: 0.9624
Epoch 10/30
11011/11011 [==============================] - 9s 799us/step - loss: 0.0440 - binary_accuracy: 0.9879 - val_loss: 0.1030 - val_binary_accuracy: 0.9619
Epoch 11/30
11011/11011 [==============================] - 10s 945us/step - loss: 0.0403 - binary_accuracy: 0.9897 - val_loss: 0.1033 - val_binary_accuracy: 0.9616
Epoch 12/30
11011/11011 [==============================] - 9s 849us/step - loss: 0.0370 - binary_accuracy: 0.9906 - val_loss: 0.1040 - val_binary_accuracy: 0.9613
Epoch 13/30
11011/11011 [==============================] - 10s 947us/step - loss: 0.0339 - binary_accuracy: 0.9922 - val_loss: 0.1039 - val_binary_accuracy: 0.9616
Epoch 14/30
11011/11011 [==============================] - 10s 884us/step - loss: 0.0313 - binary_accuracy: 0.9930 - val_loss: 0.1063 - val_binary_accuracy: 0.9597
Epoch 15/30
11011/11011 [==============================] - 11s 961us/step - loss: 0.0289 - binary_accuracy: 0.9941 - val_loss: 0.1051 - val_binary_accuracy: 0.9602
Epoch 16/30
11011/11011 [==============================] - 9s 823us/step - loss: 0.0266 - binary_accuracy: 0.9949 - val_loss: 0.1058 - val_binary_accuracy: 0.9610
Epoch 17/30
11011/11011 [==============================] - 9s 806us/step - loss: 0.0244 - binary_accuracy: 0.9959 - val_loss: 0.1075 - val_binary_accuracy: 0.9597
Epoch 18/30
11011/11011 [==============================] - 9s 817us/step - loss: 0.0227 - binary_accuracy: 0.9959 - val_loss: 0.1078 - val_binary_accuracy: 0.9616
Epoch 19/30
11011/11011 [==============================] - 9s 802us/step - loss: 0.0210 - binary_accuracy: 0.9965 - val_loss: 0.1084 - val_binary_accuracy: 0.9616
Epoch 20/30
11011/11011 [==============================] - 9s 818us/step - loss: 0.0194 - binary_accuracy: 0.9966 - val_loss: 0.1094 - val_binary_accuracy: 0.9616
Epoch 21/30
11011/11011 [==============================] - 9s 823us/step - loss: 0.0179 - binary_accuracy: 0.9971 - val_loss: 0.1109 - val_binary_accuracy: 0.9597
Epoch 22/30
11011/11011 [==============================] - 9s 803us/step - loss: 0.0166 - binary_accuracy: 0.9972 - val_loss: 0.1113 - val_binary_accuracy: 0.9616
Epoch 23/30
11011/11011 [==============================] - 9s 842us/step - loss: 0.0154 - binary_accuracy: 0.9974 - val_loss: 0.1124 - val_binary_accuracy: 0.9602
Epoch 24/30
11011/11011 [==============================] - 9s 850us/step - loss: 0.0142 - binary_accuracy: 0.9978 - val_loss: 0.1132 - val_binary_accuracy: 0.9613
Epoch 25/30
11011/11011 [==============================] - 9s 817us/step - loss: 0.0132 - binary_accuracy: 0.9983 - val_loss: 0.1144 - val_binary_accuracy: 0.9605
Epoch 26/30
11011/11011 [==============================] - 9s 848us/step - loss: 0.0122 - binary_accuracy: 0.9985 - val_loss: 0.1153 - val_binary_accuracy: 0.9605
Epoch 27/30
11011/11011 [==============================] - 10s 865us/step - loss: 0.0114 - binary_accuracy: 0.9985 - val_loss: 0.1164 - val_binary_accuracy: 0.9602
Epoch 28/30
11011/11011 [==============================] - 11s 1ms/step - loss: 0.0106 - binary_accuracy: 0.9986 - val_loss: 0.1171 - val_binary_accuracy: 0.9602
Epoch 29/30
11011/11011 [==============================] - 9s 859us/step - loss: 0.0099 - binary_accuracy: 0.9987 - val_loss: 0.1177 - val_binary_accuracy: 0.9616
Epoch 30/30
11011/11011 [==============================] - 9s 819us/step - loss: 0.0093 - binary_accuracy: 0.9988 - val_loss: 0.1188 - val_binary_accuracy: 0.9608
115/115 [==============================] - 0s 738us/step - loss: 0.1188 - binary_accuracy: 0.9608
Test accuracy: 0.9607736468315125


M_drop_fil= pd.read_csv(r"C:\Users\marti\OneDrive\AI_Lab\MCF7_DropS_Filtered_Normalised_3000_Data_train.txt", sep = ' ')


M_drop_fil.head()


M_drop_compl= pd.DataFrame(M_drop_fil.T)
M_drop_compl['label']= np.zeros(M_drop_fil.shape[1])

for x in M_drop_fil:
    if 'Normoxia' in x:
        M_drop_compl['label'][x] = 0
    if 'Hypoxia' in x:
        M_drop_compl['label'][x] = 1

M_drop_compl['label']= M_drop_compl['label'].astype(int)
M_drop_compl=  M_drop_compl.T


# call X and y
X= M_drop_fil.T
y= M_drop_compl.T['label']


X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size = 0.2, random_state=42)


mclf = RandomForestClassifier(random_state=43)
mclf.fit(X_train,y_train)

RandomForestClassifier(random_state=43)


mclf = RandomForestClassifier(random_state=43)
mclf.fit(X_train,y_train)
important_features = pd.DataFrame(mclf.feature_importances_, index = M_drop_compl.T.drop('label', axis = 1).columns).rename(columns = {0: 'Important features'})
important_features = important_features[important_features['Important features'] > 0]
important_features


X= X[important_features.index]


scores = pd.DataFrame()
scores = pd.DataFrame()
scores['accuracy_score'] = []
scores['f1_score'] = []
scores['precision_score'] = []
scores['recall_score'] = []
cross_ent = []
MSE = []


X_train, X_test, y_train, y_test = model_selection.train_test_split(X, y, test_size = 0.25, random_state=43)


nadam_model = Sequential([
    Dense(64, activation='relu', input_dim=X_train.shape[1]),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(64, activation='relu'),
    Dense(1, activation='sigmoid'),
])
nadam_model.compile(optimizer='Nadam', loss='binary_crossentropy', metrics=['binary_accuracy'])
nadam_model.fit(X_train, y_train, epochs=15, batch_size=32, validation_data=(X_test, y_test))
loss, accuracy = nadam_model.evaluate(X_test, y_test)
print(f'Test accuracy: {accuracy}')

Epoch 1/15
507/507 [==============================] - 2s 3ms/step - loss: 0.0931 - binary_accuracy: 0.9638 - val_loss: 0.0546 - val_binary_accuracy: 0.9798
Epoch 2/15
507/507 [==============================] - 1s 2ms/step - loss: 0.0293 - binary_accuracy: 0.9904 - val_loss: 0.0784 - val_binary_accuracy: 0.9737
Epoch 3/15
507/507 [==============================] - 2s 3ms/step - loss: 0.0136 - binary_accuracy: 0.9954 - val_loss: 0.0740 - val_binary_accuracy: 0.9791
Epoch 4/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0101 - binary_accuracy: 0.9963 - val_loss: 0.0778 - val_binary_accuracy: 0.9806
Epoch 5/15
507/507 [==============================] - 2s 3ms/step - loss: 0.0070 - binary_accuracy: 0.9981 - val_loss: 0.1372 - val_binary_accuracy: 0.9700
Epoch 6/15
507/507 [==============================] - 2s 3ms/step - loss: 0.0025 - binary_accuracy: 0.9992 - val_loss: 0.1171 - val_binary_accuracy: 0.9804
Epoch 7/15
507/507 [==============================] - 1s 3ms/step - loss: 6.3715e-04 - binary_accuracy: 0.9998 - val_loss: 0.1575 - val_binary_accuracy: 0.9798
Epoch 8/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0054 - binary_accuracy: 0.9982 - val_loss: 0.1064 - val_binary_accuracy: 0.9793
Epoch 9/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0044 - binary_accuracy: 0.9986 - val_loss: 0.1112 - val_binary_accuracy: 0.9778
Epoch 10/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0025 - binary_accuracy: 0.9991 - val_loss: 0.1320 - val_binary_accuracy: 0.9819
Epoch 11/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0032 - binary_accuracy: 0.9991 - val_loss: 0.1106 - val_binary_accuracy: 0.9806
Epoch 12/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0017 - binary_accuracy: 0.9996 - val_loss: 0.1286 - val_binary_accuracy: 0.9776
Epoch 13/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0014 - binary_accuracy: 0.9996 - val_loss: 0.1432 - val_binary_accuracy: 0.9795
Epoch 14/15
507/507 [==============================] - 1s 3ms/step - loss: 0.0019 - binary_accuracy: 0.9993 - val_loss: 0.1182 - val_binary_accuracy: 0.9821
Epoch 15/15
507/507 [==============================] - 1s 3ms/step - loss: 1.9680e-04 - binary_accuracy: 0.9999 - val_loss: 0.1407 - val_binary_accuracy: 0.9817
169/169 [==============================] - 0s 695us/step - loss: 0.1407 - binary_accuracy: 0.9817
Test accuracy: 0.9816904067993164

	"output.STAR.PCRPlate1A10_Normoxia_S123_Aligned.sortedByCoord.out.bam"	"output.STAR.PCRPlate1A4_Hypoxia_S8_Aligned.sortedByCoord.out.bam"	"output.STAR.PCRPlate1A5_Hypoxia_S108_Aligned.sortedByCoord.out.bam"	"output.STAR.PCRPlate1A7_Normoxia_S113_Aligned.sortedByCoord.out.bam"	...	"output.STAR.PCRPlate4H10_Normoxia_S210_Aligned.sortedByCoord.out.bam"	"output.STAR.PCRPlate4H2_Hypoxia_S199_Aligned.sortedByCoord.out.bam"
"WASH7P"	0	0	0	0	...	0	0
"CICP27"	0	0	0	0	...	0	0
"DDX11L17"	0	1	0	0	...	0	0
"WASH9P"	0	0	0	1	...	1	1
"OR4F29"	2	0	1	0	...	0	0

	S123	S26	S97	S104	S4	S8	S108	S11	S113	S119	...	S243	S193	S198	S232	S204	S210	S214	S199	S205	S236
count	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	...	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000	23396.000000
mean	99.565695	207.678278	9.694734	150.689007	35.700504	47.088434	152.799453	135.869422	38.363908	45.512139	...	76.361771	105.566593	54.026116	29.763806	28.905411	104.740725	35.181569	108.197940	37.279962	76.303855
std	529.532443	981.107905	65.546050	976.936548	205.885369	545.367706	864.974182	870.729740	265.062493	366.704721	...	346.659348	536.881574	344.068304	186.721266	135.474736	444.773045	170.872090	589.082268	181.398951	369.090274
min	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
50%	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	...	1.000000	0.000000	0.000000	0.000000	0.000000	2.000000	0.000000	0.000000	0.000000	1.000000
75%	51.000000	125.000000	5.000000	40.000000	22.000000	17.000000	81.000000	76.000000	22.000000	18.000000	...	56.000000	67.000000	29.000000	18.000000	19.000000	76.000000	24.000000	68.000000	22.000000	44.000000
max	35477.000000	69068.000000	6351.000000	70206.000000	17326.000000	47442.000000	43081.000000	62813.000000	30240.000000	35450.000000	...	19629.000000	30987.000000	21894.000000	13457.000000	11488.000000	33462.000000	15403.000000	34478.000000	10921.000000	28532.000000

	mean	std	min	max
WASH7P	0.040609	0.316833	0.0	3.0
CICP27	0.086294	0.481987	0.0	5.0
DDX11L17	0.380711	1.116771	0.0	8.0
WASH9P	0.187817	0.580512	0.0	4.0
OR4F29	0.091371	0.352508	0.0	2.0
...	...	...	...	...
MT-TE	14.472081	19.036096	0.0	212.0
MT-CYB	2023.000000	1375.704899	83.0	6760.0
MT-TT	19.959391	21.538489	0.0	154.0
MT-TP	42.020305	35.000723	0.0	267.0
MAFIP	3.345178	4.108557	0.0	24.0

	P-value
Gene
CA9	1.548012e-32
ANGPTL4	5.527006e-32
NDRG1	6.695692e-32
EGLN3	8.133323e-32
BNIP3P1	2.672800e-26
KCTD11	3.891886e-25
LOXL2	4.807956e-25
ASB2	6.391086e-24
ENO2	3.995859e-23
PPP1R3G	1.702522e-20
CCNG2	7.751523e-19
DANT1	1.265791e-18
PNRC1	1.183295e-16
BARX1	1.293553e-16
HK2	2.989836e-16

	DDX11L17	WASH9P	OR4F29	MTCO2P12	MTATP8P1	LINC01409	FAM87B	LINC01128	LINC00115	KLHL17	...	C21orf58	DIP2A-IT1	DSTNP1	MT-TF	MT-TI	MT-TQ	MT-TM	MT-TD	MT-TR	MAFIP
0
1	-0.620925	-0.620925	-0.620925	-0.620925	-0.620925	-0.620925	-0.620925	-0.620925	-0.620925	2.423598	...	2.087125	-0.620925	-0.620925	-0.620925	0.072222	0.072222	0.988513	0.072222	-0.620925	-0.620925
1	0.138874	-0.554273	-0.554273	-0.554273	-0.554273	-0.554273	-0.554273	-0.554273	-0.554273	-0.554273	...	-0.554273	-0.554273	-0.554273	-0.554273	0.544339	0.544339	-0.554273	-0.554273	-0.554273	1.748312
1	-0.921032	-0.921032	-0.921032	0.870727	-0.921032	-0.921032	-0.921032	-0.921032	-0.921032	-0.921032	...	2.689886	-0.921032	1.276192	-0.227885	-0.227885	0.465262	0.465262	-0.227885	0.177580	0.688406
1	1.679592	-0.517633	0.580980	-0.517633	-0.517633	-0.517633	-0.517633	2.978875	-0.517633	-0.517633	...	2.740464	-0.517633	-0.517633	-0.517633	-0.517633	0.175514	0.868662	-0.517633	-0.517633	-0.517633
1	0.635405	-0.463208	0.229940	-0.463208	-0.463208	-0.463208	-0.463208	2.714846	-0.463208	-0.463208	...	-0.463208	0.229940	-0.463208	-0.463208	-0.463208	-0.463208	0.229940	0.229940	-0.463208	-0.463208

	Cell Line	PCR Plate	Pos	Condition	Hours	Cell name	PreprocessingTag	ProcessingComments
Filename
output.STAR.PCRPlate1A10_Normoxia_S123_Aligned.sortedByCoord.out.bam	HCC1806	1	A10	Normo	24	S123	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
output.STAR.PCRPlate1A12_Normoxia_S26_Aligned.sortedByCoord.out.bam	HCC1806	1	A12	Normo	24	S26	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
output.STAR.PCRPlate1A1_Hypoxia_S97_Aligned.sortedByCoord.out.bam	HCC1806	1	A1	Hypo	24	S97	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
output.STAR.PCRPlate1A2_Hypoxia_S104_Aligned.sortedByCoord.out.bam	HCC1806	1	A2	Hypo	24	S104	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
output.STAR.PCRPlate1A3_Hypoxia_S4_Aligned.sortedByCoord.out.bam	HCC1806	1	A3	Hypo	24	S4	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts

	Cell Line	PCR Plate	Condition	Hours
Cell name
S123	HCC1806	1	Normo	24
S26	HCC1806	1	Normo	24
S97	HCC1806	1	Hypo	24
S104	HCC1806	1	Hypo	24
S4	HCC1806	1	Hypo	24
...	...	...	...	...
S210	HCC1806	4	Normo	24
S214	HCC1806	4	Normo	24
S199	HCC1806	4	Hypo	24
S205	HCC1806	4	Normo	24
S236	HCC1806	4	Normo	24

	S123	S26	S97	S104	S4	S8	S108	S11	S113	S119	...	S243	S193	S198	S232	S204	S210	S214	S199	S205	S236
MTCO1P12	1740	2340	22	1020	303	587	2572	1004	443	779	...	639	790	644	386	297	961	372	1128	258	980
MTATP6P1	974	1642	8	653	91	161	1242	438	171	507	...	422	496	253	171	213	642	175	638	162	467
NOC2L	285	474	77	1987	162	0	955	648	298	9	...	581	860	156	205	360	1352	166	558	522	286
ISG15	256	897	17	919	126	247	267	530	266	121	...	146	346	180	72	42	1142	40	394	58	275
SDF4	155	376	51	436	112	311	89	426	162	112	...	375	535	100	140	95	496	41	439	328	134
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
MT-ND4L	6060	4281	82	1076	279	561	4615	1208	406	1418	...	1450	1393	770	495	442	1977	441	2171	451	1621
MT-ND4	17732	15624	200	6588	1028	1768	14891	4139	1738	5964	...	4121	3447	2571	1779	1395	6273	1446	7026	1201	5341
MT-ND5	3852	7457	93	1479	303	650	6533	1271	589	1574	...	1362	1993	884	746	555	2926	688	2999	611	2460
MT-ND6	900	1439	12	234	33	79	1281	250	65	331	...	310	344	163	119	53	423	114	486	75	513
MT-CYB	4208	6491	25	4819	310	695	2885	1552	366	1829	...	1119	1429	808	999	916	3719	984	2256	981	2039

	0	1	2	3	4	5	6	7	8	9	...	165	166	167	168	169	170	171	172	173	174
true_label
1	0.962692	5.231543	-14.155139	3.840035	-14.940994	-14.564328	14.718075	-7.491317	5.269369	-1.142132	...	-5.961566	-0.392383	5.200690	-6.646728	-10.589341	-2.859760	-6.064934	-2.592688	-0.994372	5.323559
1	52.151008	31.430537	4.498231	4.049625	1.064632	-5.853009	21.913842	-6.861485	-13.361533	-1.090423	...	0.324911	0.506792	1.358075	2.348563	0.318519	1.279852	5.125902	1.990882	-0.891389	2.583927
1	-9.738035	12.012388	-1.495742	0.775028	-29.889254	22.762295	8.425296	8.607963	7.096413	20.018429	...	-0.510999	2.120077	-2.478407	0.054196	0.003741	2.024505	-0.824480	1.357065	0.764387	0.955908
1	10.862382	10.430442	-36.099757	8.437602	-16.332751	-12.874535	14.357779	0.611470	3.099180	-8.453487	...	-0.536031	1.996279	-0.072826	-3.186667	1.458440	-3.180844	-0.840788	-1.354318	0.136751	0.106143
1	15.290039	11.436769	-28.681119	-9.996548	-1.280457	-8.713127	5.457373	-3.122674	4.393358	-3.889854	...	-7.938023	5.364352	1.731742	3.018310	3.641144	-1.349306	-6.644473	-5.412445	9.468888	0.551855
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
0	2.247701	-10.188355	-2.473914	-10.871130	-9.772319	-13.634096	6.703722	-9.991245	-4.600410	-0.131346	...	8.561964	4.613611	-6.630280	-7.304093	-19.172097	-13.929491	-17.249690	7.889787	-5.370012	-4.802503
0	-33.687925	-8.954735	22.868723	-6.982392	1.937430	4.240745	10.789240	-2.336619	7.188888	8.124489	...	1.731812	0.820602	1.810733	0.788486	-1.370114	1.556318	-2.274760	-0.334841	0.984749	-2.139395
0	-13.397524	-7.061122	-10.881664	-17.257769	13.571037	-3.832316	9.535070	-2.044857	-0.723982	-7.276656	...	18.008331	16.213082	3.697997	6.158550	7.452212	-1.607394	13.498858	0.800144	3.415102	2.525905
0	5.456924	-9.911785	-5.662801	-14.821987	6.327714	-3.165563	3.666687	-5.905843	8.005732	0.072053	...	13.233844	-11.303779	2.515552	-9.762822	-10.406417	16.582153	-18.012140	-12.585750	-19.306883	10.312608
0	-8.897698	-27.228986	5.126364	4.093408	-12.922448	0.815748	-5.653270	9.842139	-0.926639	0.614727	...	7.357868	8.632228	2.375766	7.846683	-4.405730	2.691538	0.812541	1.125435	-0.284296	-2.262204

	Important features
KLHL17	0.000374
RN7SL657P	0.000169
TMEM52	0.000088
HES5	0.001027
TNFRSF14-AS1	0.000079
...	...
DYRK1A	0.000229
FAM3B	0.000083
GATD3A	0.000611
MT-TF	0.001565
MAFIP	0.000816

	Important features
CA9	0.029320
LOXL2	0.024512
PPP1R3G	0.021394
KCTD11	0.017708
ANGPTL4	0.015491
NDRG1	0.013857
ENO2	0.011364
EGLN3	0.010957
PAG1	0.010453
BNIP3P1	0.010273
PFKFB4	0.010032
DANT1	0.009983
CCNG2	0.009843
ASB2	0.009556
GPR146	0.008232

	Important features
DDIT4	0.024601
ANGPTL4	0.016027
CALML5	0.000184
CCNB1	0.000275
IGFBP3	0.000265
...	...
APLF	0.000271
ST14	0.000409
INSYN2B	0.000273
NMB	0.000241
NEAT1	0.001262

	S57	S58	S59	S60	S79	S81	S130	S131	S103	S104	...	S382	S383	S355	S357	S358	S359	S360	S379	S380	S381
CYP1B1	343	131	452	27	5817	3841	9263	21543	1013	53	...	7890	4512	160	351	327	196	504	34565	20024	5953
CYP1B1-AS1	140	59	203	7	2669	1565	3866	9113	459	22	...	3647	2035	75	138	130	102	238	13717	7835	2367
CYP1A1	0	0	0	0	0	79	238	443	0	0	...	86	1654	0	0	0	1	0	11274	563	522
NDRG1	0	1	0	0	654	1263	2634	540	0	13	...	481	1052	0	0	54	243	62	1263	925	1572
DDIT4	386	289	0	288	2484	2596	1323	2044	36	204	...	3692	2410	800	1	189	266	417	4256	12733	2275

	Cell Line	Lane	Pos	Condition	Hours	Cell name	PreprocessingTag	ProcessingComments
Filename
S28	MCF7	output.STAR.1	A10	Hypo	72	S28	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
S29	MCF7	output.STAR.1	A11	Hypo	72	S29	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
S30	MCF7	output.STAR.1	A12	Hypo	72	S30	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
S1	MCF7	output.STAR.1	A1	Norm	72	S1	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts
S2	MCF7	output.STAR.1	A2	Norm	72	S2	Aligned.sortedByCoord.out.bam	STAR,FeatureCounts

	Important features
CYP1B1	0.008707
CYP1B1-AS1	0.023430
CYP1A1	0.000978
NDRG1	0.022070
DDIT4	0.007328
...	...
LINC02861	0.000465
SGMS1-AS1	0.000191
ABAT	0.000874
FZD1	0.000448
PLEKHN1	0.000713

	KLHL17	RN7SL657P	TMEM52	HES5	TNFRSF14-AS1	TP73-AS1	DFFB	AJAP1	NPHP4	KCNAB2	...	BAGE2	CYYR1	GPX1P2	SIM2	HLCS	DYRK1A	FAM3B	GATD3A	MT-TF	MAFIP
0
1	2.423598	-0.620925	-0.620925	-0.620925	-0.620925	2.746371	-0.620925	-0.620925	4.207389	0.477687	...	-0.620925	-0.620925	-0.620925	1.458517	3.473420	0.988513	-0.620925	-0.620925	-0.620925	-0.620925
1	-0.554273	-0.554273	-0.554273	-0.554273	2.664603	-0.554273	-0.554273	-0.554273	-0.554273	1.930633	...	-0.554273	-0.554273	-0.554273	-0.554273	3.694222	1.237486	-0.554273	-0.554273	-0.554273	1.748312
1	-0.921032	-0.921032	-0.921032	-0.921032	-0.921032	1.158409	3.189842	3.998949	2.170010	-0.921032	...	0.465262	-0.921032	-0.921032	-0.921032	3.632845	1.476863	2.716554	-0.921032	-0.227885	0.688406
1	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	0.580980	-0.517633	...	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633	-0.517633
1	-0.463208	-0.463208	-0.463208	-0.463208	2.101742	-0.463208	-0.463208	-0.463208	-0.463208	-0.463208	...	0.229940	-0.463208	-0.463208	2.627835	-0.463208	-0.463208	-0.463208	-0.463208	-0.463208	-0.463208
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
0	-0.560286	-0.560286	-0.560286	-0.560286	-0.560286	-0.560286	-0.560286	-0.560286	1.519156	-0.560286	...	-0.560286	-0.560286	-0.560286	2.004664	-0.560286	-0.560286	-0.560286	-0.560286	-0.560286	0.132862
0	3.586641	-1.192482	-1.192482	0.599277	0.753428	0.193812	-1.192482	3.870113	2.208715	2.777810	...	0.193812	-1.192482	-1.192482	-1.192482	-1.192482	1.372467	-1.192482	-1.192482	1.004742	0.193812
0	-0.745941	-0.745941	1.451284	-0.745941	-0.745941	-0.745941	-0.745941	1.651954	-0.745941	-0.745941	...	-0.052794	-0.745941	-0.745941	2.472935	2.809407	1.333501	1.819008	-0.745941	-0.745941	-0.745941
0	-0.579741	-0.579741	-0.579741	-0.579741	-0.579741	-0.579741	-0.579741	1.722844	-0.579741	-0.579741	...	-0.579741	-0.579741	-0.579741	-0.579741	2.975607	0.806553	-0.579741	-0.579741	0.113406	0.113406
0	-0.961269	-0.961269	2.174225	-0.961269	0.648169	-0.961269	2.727610	0.648169	1.436626	1.983170	...	-0.268122	-0.961269	-0.961269	-0.961269	0.425025	-0.961269	-0.961269	-0.961269	-0.268122	0.830490

	AAAAAACCCGGC_Normoxia	AAAACCGGATGC_Normoxia	AAAACGAGCTAG_Normoxia	AAAACTTCCCCG_Normoxia	AAAAGCCTACCC_Normoxia	AAACACAAATCT_Normoxia	AAACCAAGCCCA_Normoxia	AAACCATGCACT_Normoxia	AAACCTCCGGCT_Normoxia	AAACGCCGGTCC_Normoxia	...	TTTTCTGATGGT_Hypoxia	TTTTGATTCAGA_Hypoxia	TTTTGCAACTGA_Hypoxia	TTTTGCCGGGCC_Hypoxia	TTTTGTTAGCCT_Hypoxia	TTTTTCCGTGCA_Hypoxia	TTTTTGCCTGGG_Hypoxia	TTTTTGTAACAG_Hypoxia	TTTTTTTGAATC_Hypoxia
H1-5	2	2	5	1	0	0	0	0	1	0	...	0	1	0	2	1	0	0	3	1
MALAT1	3	3	2	3	12	3	1	2	0	0	...	3	1	1	1	4	4	1	3	6
MT-RNR2	0	0	0	0	0	0	0	0	0	1	...	1	2	2	2	0	1	0	1	0
ARVCF	0	0	0	0	0	0	0	0	0	0	...	0	0	0	0	0	0	0	0	0
BCYRN1	0	1	1	0	0	1	1	2	0	3	...	1	1	0	1	1	0	1	0	0

	AAAAAACCCGGC_Normoxia	AAAACCGGATGC_Normoxia	AAAACGAGCTAG_Normoxia	AAAACTTCCCCG_Normoxia	AAAAGCCTACCC_Normoxia	AAACACAAATCT_Normoxia	AAACCAAGCCCA_Normoxia	AAACCATGCACT_Normoxia	AAACCTCCGGCT_Normoxia	AAACGCCGGTCC_Normoxia	...	TTTTCTGATGGT_Hypoxia	TTTTGATTCAGA_Hypoxia	TTTTGCAACTGA_Hypoxia	TTTTGCCGGGCC_Hypoxia	TTTTGTTAGCCT_Hypoxia	TTTTTACCAATC_Hypoxia	TTTTTCCGTGCA_Hypoxia	TTTTTGCCTGGG_Hypoxia	TTTTTGTAACAG_Hypoxia	TTTTTTTGAATC_Hypoxia
count	3000.00000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.00000	3000.000000	3000.000000	...	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000	3000.000000
mean	0.02900	0.041667	0.024333	0.021667	0.029667	0.020000	0.036000	0.02600	0.034000	0.029333	...	0.043000	0.049667	0.037000	0.047667	0.057000	0.023333	0.041667	0.041667	0.043333	0.040000
std	0.23276	0.309778	0.231860	0.189409	0.323761	0.170126	0.250449	0.23525	0.231362	0.218683	...	0.271739	0.319219	0.279864	0.259648	0.304053	0.214797	0.236536	0.285116	0.267356	0.282418
min	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
25%	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
50%	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
75%	0.00000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.00000	0.000000	0.000000	...	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000	0.000000
max	5.00000	9.000000	7.000000	4.000000	12.000000	3.000000	4.000000	6.00000	4.000000	4.000000	...	4.000000	7.000000	7.000000	4.000000	5.000000	4.000000	4.000000	5.000000	5.000000	6.000000

	Important features
H1-5	0.002466
MALAT1	0.009115
MT-RNR2	0.003885
ARVCF	0.000025
BCYRN1	0.015336
...	...
SCCPDH	0.000099
NTAN1	0.000096
CLIP2	0.000051
DUSP23	0.000155
ZNF682	0.000018

	Important features
MALAT1	0.007348
MT-RNR2	0.033106
NEAT1	0.001162
H1-5	0.001020
TFF1	0.037403
...	...
BRWD1-AS2	0.000008
RPS19BP1	0.000180
AUNIP	0.000024
TNK2	0.000050
SUDS3	0.000394

AI Lab project

Introduction¶

Cancer and Genomic research¶

Hypoxia and its role in cancer recognition¶

Target of this research¶

Preamble¶

Overview over data¶

Steps of our analysis¶

Imports¶

Analysis of HCC1806_SmartS_Unfiltered¶

Data Exploration¶

UNSUPERVISED ANALYSIS¶

Data Cleaning and Quality Control¶

Cell Filtering¶

Step 1: Outliers detection¶

Step 2: Mitochondrial reads¶

Gene Filtering¶

First approach: Duplicates¶

Second Approach: Coefficient of variation¶

Third Approach: Sparsity¶

Normalisation and Scaling¶

First Normalization Approach¶

Second Approach¶

MODEL FOR OUR DATA¶

UNSUPERVISED LEARNING ANALYSIS¶

SUPERVISED LEARNING¶

On the whole dataset¶

RANDOM FOREST¶

FEATURE SELECTION¶

SVM¶

On the reduced dataset with RandomForest¶

On the reduced dataset with PCA¶

CHOSEN MODEL¶

From now on...¶

HCC1806_SmartS_Filtered_Normalised¶

About EDA¶

ANALYSIS¶

Unsupervised Learning¶

K-means clustering¶

Feature selection¶

Model Selection¶

MCF7_SmartS_Filtered_Normalised¶

SOME VISUALIZATION¶

Data Preparation¶

FEATURE SELECTION¶

Random forest¶

Logistic Regression¶

Model with important features data¶

HCC1806_DropS_Filtered_Normalised¶

Exploratory Data Analysis(EDA)¶

More details about the distributions:¶

Correlation¶

Target Variables¶

Visualization of the data¶

Supervised Learning Model¶

Feature Selection¶

Neural Network¶

MCF7_DropS_Filtered_Normalised¶

Supervised Learning Model¶

Feature selection¶

Neural Network¶

CONCLUSIONS¶

BIBLIOGRAPHY¶